{'fraction_non_alphanumeric': 0.04151983725557767, 'fraction_numerical': 0.01529619707193437, 'mean_word_length': 4.084958453450907, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dual Subtractions

We propose a novel local subtraction scheme for the computation of Next-to-Leading Order contributions to theoretical predictions for scattering processes in perturbative Quantum Field Theory. With respect to well known schemes proposed since many years that build upon the analysis of the real radiation matrix elements, our construction starts from the loop diagrams and exploit their loop-tree dual representation. Our scheme implements exact phase space factorization, handles final state as well as initial state singularities and is suitable for both massless and massive particles.


Introduction
The success of the physics studies at the LHC, culminated with the discovery of the Higgs boson [1,2] at CERN [3,4], has its roots in the deep level of understanding reached in both experimental and theoretical aspects of the physics of hadronic collisions. On the theory side, the computation of scattering amplitudes including higher order perturbative corrections plays the main role. The calculation of tree-level and one-loop matrix elements is nowadays fully automated in very efficient ways. Furthermore, well known schemes for the treatment of their infrared and collinear behaviour have been formulated since decades, like the Catani-Seymour and the Frixione-Kunszt-Signer local schemes [5][6][7]. For these reasons the Next to Leading Order (NLO) computation is considered a completely solved problem. On the other hand, new precise exclusive measurements will be available soon, as the ongoing monumental campaign of data collection goes on, so that, for a meaningful comparison among theory and experiment, Next-to-Next-to Leading Order (NNLO) calculations are required. At the NNLO, one is faced with the construction of appropriate integration schemes that allows for the cancellation of non integrable infrared and collinear singularities across contributions that live on three different phase spaces. Several techniques have been developed to address this problem, such as q T [8] and N -jettiness [9,10] slicing methods, subtraction schemes like antenna [11][12][13][14][15][16][17][18][19], CoLoRFulNNLO [20][21][22][23][24][25][26][27][28][29], residue-improved [30][31][32][33], nested soft-collinear [34][35][36][37][38] and projection-to-Born [39]. Also, other approaches are under development [40][41][42][43]. Although there is still not the same level of automation as for the NLO, progress on the subtraction of the singularities is going quite fast.
However, the bottleneck for NNLO computations is represented by the calculation of multi-loop amplitudes with several massive internal and external particles. This field of research is very active, although it is still not clear if the procedures proposed so far can really simplify the job. The most efficient paradigm followed for decades keeps separated the two problems mentioned above. Reduction of tensor integrals to a base of master integrals and computation of the master integrals from one side, and the construction of a subtraction scheme based on the analysis of the radiation matrix elements, on the other.
An interesting new path was proposed in [44,45]. The computation of loop diagrams is turned into phase space integrals thanks to a Loop-Tree duality (LTD) theorem proved in the same papers. The possibility to follow such a strategy is of course very interesting because the high mathematical complexity accompanying the analytic computation of higher loop calculation and that of building a subtraction scheme, might simply not be there, being replaced by the numerical integration of properly defined integrands. Such integrations are of course non trivial, but one can think that computer science has already developed a large set of tools to address the associated technicalities. This path has been pursued in [46][47][48] at one-loop and in [49] at two loop as examples. Another point of view offered by the LTD theorem is the possibility to simplify the direct numerical computation of multi-loop amplitudes. Non-trivial numerical applications of LTD have been performed in [50][51][52]. Furthermore, in [53][54][55] the better numerical behaviour of LTD due to causality has been conjectured and demonstrated. Further steps in this direction have been done in [56][57][58]. In particular, in [56] a general strategy is derived for the subtraction of the divergences by one and two loop QED amplitudes with the effect that the finite remnants can be numerically evaluated with relatively small effort.
In the present paper, we will consider a slightly different perspective from the ones outlined above. We investigate the possibility to extract the divergences from the loop amplitude in a way that builds subtractions for the real radiation contribution. This path is somewhat opposite to the usual strategy to build subtractions, that is based on the analysis of the divergences in the real sector, nevertheless universality and cancellation of course must hold in both directions. The opportunity to follow such a path is offered by LTD, but we will not apply it directly to one-loop Feynman diagrams. We will first consider the reduction of the tensor integrals to scalar integrals and then apply LTD to a small class of scalar integrals, made by namely just two integrals, one triangle and one massless bubble (in massless QCD). Note that, to build a proper subtraction we will also need an appropriate integrand expression for the wave-function renormalization. We recall the available formulas for such wave-function contributions for the case of fermions and compute the one for the gluon. We note that a procedure to build an expression for the gluon wave-function renormalization constant at the integrand-level has also been already presented in [59]. In that work, counterterms for the virtual contribution are built and shown to cancel the singularities of the counterterms for the real part. On the other hand, our approach consists in building counterterms for the real contributions through their direct extraction from the virtual part.
The paper is organized as follows: in Section 2 we discuss the singular behaviour of oneloop matrix elements while in Section 3 we report the basics of LTD and use them to extract divergent counterterms for the real radiation from final state massless particles. Then, in Sections 4 and 5 we follow the same path to build subtractions for the real contribution in the cases of radiation from massive final state particles and from initial state particles, respectively. In Section 6 we show a small collection of applications and in Section 7 we comment our construction analyzing some consequences and give our conclusion. Finally, we include two appendices, one with about the counting of the dual counterterms (Appendix A) and a second one (Appendix B) collecting all the formulae for the numerical implementation of the dual subtraction scheme for NLO QCD corrections to scattering processes involving any number of massless external partons and non coloured particles.

Singular behavior of one-loop matrix elements
We start by considering the easiest situation that is the one of processes with a colourless initial state and m massless partons in the final state. Although most of the discussion in this section already applies to the case of final states with massive particles, and that of coloured particles in the initial state, in the present section we limit the discussion to the simpler case mentioned above and postpone the extension to the other cases to dedicated following sections. Our first step is to consider the full reduction of a one-loop amplitude to scalar integrals. As it is well known a reduced amplitude has the following form in terms of scalar boxes, triangles, bubbles and tadpoles where the indices i, j, k, l run over all possible external momenta {p} and their sums. It is worth to note that up to this integral transformation from its expression in terms of Feynman diagrams, the amplitude has retained exactly the same meaning and values at all orders in the dimensional parameter that we also keep in both the scalar integrals and their coefficients. The same is obviously true also if we use different sets or combinations of the scalar integrals. A particularly convenient transformation consists in expressing the scalar boxes in terms of its six dimensional version (D 6−2ε 0 ) or, equivalently, in terms of the rank two form factor that multiplies the metric tensor in a covariant decomposition (that we conventionally denote by D 00 ). After this exchange, Eq.(2.1) becomes The D 00 function consists of a linear combination of the relative box function and the four triangle functions that are obtained pinching one denominator at the time. This is the reason why we have not changed the coefficients of the scalar bubble and tadpole integrals. Now we analyze the infrared and collinear behaviour of the reduced amplitude. First, we remind that, irrespective from the value of both the internal and external masses, the D 00 function is completely finite with respect to the virtual loop integration, being in particular free of both Ultra-Violet (UV) and Infra-Red and collinear (IR) poles 1 .
Moving to triangle functions C ijk 0 , there are two kinds of them which have IR divergences. The first one, shown in Fig.(1a), has an outgoing momentum p i connected, through a massless propagator, to a combination k of at least two other external momenta; we can call it generically C i (k). The second type is shown in Fig.(1b) and has two outgoing momenta, p i and p j , connected by a massless propagator; for massless partons, it corresponds to the scalar three-point function C 0 (p i , p j ). Up to this point, the analysis equally applies also to the case of massive external partons. The real part of the triangular functions is given, in the massless case, by .
(2.4) 1 This property can be deduced by combining the scalar integrals of the reduced D00 function or by inspecting its direct calculation in terms of Feynman parameters. An alternative proof based on the analysis of the soft and collinear behaviour of rank two tensor box integrals can be found in [60]. The red colour indicates a massless propagator while the blue and the green colour going from an internal to an external line means that the mass of the propagator is the same of the external (on-shell) particle with the same colour.
Finally, massless scalar two-point function B ij 0 depending on a single massless external momentum p i , bring IR divergences while we assume that massless tadpoles give no contributions to IR poles (i.e. we take them as vanishing). The structure of the singularities of one-loop amplitudes is completed by adding the ones brought by the renormalization procedure. While coupling renormalization involves only UV poles, wave-function renormalization of the external particles introduces both UV and IR poles. These contributions can only bring single poles and are proportional to the whole leading order matrix element. Now that we have analyzed the IR poles of Eq.(2.2), we can compare the results with the general formula for the singular behaviour of the interference among the renormalized one-loop amplitude A (1,R) m and the tree-level one A (0) m [61,62], that is (2.6) where the usual definition of color charge operators and color correlated matrix elements is implied and the constants γ i are given by By inspecting Eq.(2.6) we observe that there are no single poles with a logarithm of an invariant formed with more than two external momenta (2 p i · k). These poles are carried on individually and exclusively by the triangle integrals C i (k). Thus we deduce that for full amplitudes the coefficients of the C i (k) scalar integrals in Eq.(2.2) must vanish 2 . It is then clear that the double pole and the logarithmic part of the single pole in Eq. (2.6) are produced by the triangle functions C 0 (p i , p j ). Furthermore, the coefficients of these functions (there is one for every pair of the external partons as for the terms of the double sum in Eq.(2.6)) obtained summing over all possible Feynman diagrams contributing to a specific full one-loop process is while the remaining single poles (the non-logarithmic ones) come from the bubbles B 0 (p i ) and the wave-function renormalization counterterms ∆Z(p i ). The latters are necessary to obtain a fully local cancellation of IR singularities between the real cross section and the dual cross section that we are going to build. The knowledge of the IR poles coming from ∆Z(p i ) and Eq.(2.6) allow to determine the coefficients of the bubble integrals B 0 (p i ), that turn out to be where C q = Cq = C F and C g = C A . The triangles C 0 (p i , p j ), the bubbles B 0 (p i ) and the renormalization counterterms ∆Z(p i ) represent the key ingredients for the construction of the dual subtraction scheme. In fact, our prescription to find all the necessary dual counterterms imposes to: • collect triangles, bubbles and wave-function renormalization counterterms; • construct their dual representation by application of LTD; • extract their IR singular part by properly selecting the loop integration domains.
Following these steps, in the next section we build a set of dual subtractions.

Loop-Tree duality and counterterms
In this Section we show how to extract the singularities from the virtual amplitude constructing counterterms. We first recall the basics of LTD in section 3.1. In particular, we list the main results that we need for our construction and refer the interested reader to the seminal paper [44] where they have been derived for the first time, and also to the useful applications in [46,48], from which we borrow most of the discussion. Then, in section 3.2 we will make a short summary of the relevant well know limiting behaviours of tree-level amplitudes. In the subsequent section 3.3 we build the dual subtractions and we prove the local cancellation of the singularities in the real matrix element in section 3.4. In section 3.5 we present the results of the integrated subtractions and show the agreement with the general structure of the singularities expressed in Eq.(2.6).

Loop-Tree duality
Let us consider a generic N -particle scalar loop, which is shown in figure 2. The external momenta are denoted as p 1 , . . . , p N and are taken as outgoing, q is the counter-clockwise flowing loop momentum, and q i are the momenta of the internal lines. Momentum conservation reads while the internal momenta q i are related to q and the external momenta by which implies the choice q N = q. With this notation, the N -particle scalar one-loop integral is given by 3 are the Feynman propagators and we have used the shorthand notation µ being an arbitrary energy scale used to restore the physical dimensions of the integral. The LTD theorem states that the loop integral in Eq.(3.3) has the following dual representation (3.5) 2 We remind here that the Ci(k) triangle functions are still contained in the reduction of the D00 functions, and that their poles will anyway cancel with those of the other scalar integrals in that reduction. 3 Adding a numerator to Eq.(3.3) one can extend the discussion to any kind of integral.
whereδ(q i ) ≡ 2πiθ(q 0 i )δ(q 2 i − m 2 i ) sets the internal momentum q i on-shell and are dual propagators. Here η is an arbitrary time-like or light-like vector with positive definite energy which can be chosen as η ≡ (1, 0). The only difference between the dual propagator and the Feynman propagator G F lies in the different i0 prescription which regularizes the singularity of the right-hand side of Eq.(3.6). Using η ≡ (1, 0), we see that this dual i0 prescription depends on the sign of q 0 j − q 0 i which in turn is related, through Eq. (3.2), to the energy components of the external momenta To the right-hand side of Eq.(3.5) we find N different integrals. In each of them, one Feynman propagator is removed, the corresponding internal momentum is set on-shell and the remaining Feynman propagators are turned into dual propagators. The role of the on-shell internal momentum is passed through the N internal momenta q 1 , . . . , q N without repetition, so that each momentum q i is on-shell in one and only one integral. Moreover, in each integral we can change the integration variable in order to integrate over the onshell momentum itself. If q i is the on-shell momentum, this is done by mean of a simple translation suggested by Eq. (3.2). This means that the dual representation in Eq.(3.5) allows us to replace a one-loop integral over an off-shell virtual momentum by a linear combination of integrals over on-shell virtual momenta. Let us now consider the dual representation of a virtual amplitude, and suppose we have built a dual subtraction from it. In the framework of the subtraction method, we should integrate the dual cross section together with the real one. Since they live on different phase spaces, we need a mapping (Φ m , q i ) ↔ Φ m+1 connecting the virtual to the real phase space. To be precise, we will use different mappings for different dual counterterms. This means that, in a Monte Carlo integration, we can generate a point in the real phase space and then obtain, through the application of each mapping, as many virtual configurations as the number of counterterms. Moreover, in order to achieve a local cancellation of divergences, we require that IR singularities of the dual integrations are mapped into soft and collinear singularities of the (m + 1)-particle matrix element squared. Considering the emitter-spectator pair (i, j), we always denote the internal momenta of three-point function C 0 (p i , p j ) as in Fig.(3a), namely q 1 is the virtual momentum connecting the external particles, while q 2 = q 1 + p j and q 3 = q 1 − p i . The dual contributions where q 1 is set on-shell by LTD bring IR singularities that can be associated with a real sector parton p c being emitted from a parton of momentum p a and absorbed by a spectator p b . To this end, we make use of the following momentum mapping [46] where α (1) abc ≡ s ac /s abc , s ac = (p a + p c ) 2 (that in the massless case reduces to s ac = 2p a · p c ), s abc = (p a + p b + p c ) 2 , q 1 is the virtual on-shell loop momentum, and all the other momenta are unchanged. The definition domain of Eq.(3.9) is chosen to be the region where s ac < s bc . This selects a compact subset of the loop momentum phase space, which will be parametrized later. The momentum mapping in Eq.(3.9) automatically verifies p 1 + p 2 + p 3 = p 1 + p 2 (momentum conservation) and p 2 1 = p 2 2 = 0 (on-shell conditions). It is useful to consider also the inverse momentum mapping, which is given by and is the one explicitly used in [46]. The soft emission p c → 0 is mapped into q 1 → 0, while the collinear configuration p a ||p c is given by q 1 ||p i , as can be deduced by noticing that α (1) ij → 0 when q 1 · p i → 0. Indeed one has that Also, Eq.(3.9) clearly shows that, when p a and p c become collinear, the momentum p i is mapped into the emitter momentum p a + p c . In order to properly relate the integration over Φ m+1 to the one over (Φ m , q i ), we need to compute the jacobian of the transformation in Eq. (3.9). To do this, let us move to the center-of-mass frame of p i and p j , where the loop momenta q k can be assigned in a spherical coordinate system whose zenith is the direction of p i Here, ξ k √ s ij /2 is the energy of q k , v k ∈ [0, 1] is related to the cosine of the polar angle by cos θ k ≡ 1 − 2v k and ϕ k is the azimuthal angle. We can now express the momentum mapping in Eq.(3.9) as a relation between the dimensionless variables (ξ 1 , v 1 ) and the kinematic invariants y mn ≡ s mn /s abc , obtaining (3.13) The collinear limit y ac → 0 is mapped into v 1 → 0, while the soft one is approached when ξ 1 → 1. The jacobian related to Eq.(3.13) is given by . (3.14) Also, by using the inverse of Eq.(3.13) we can find the boundaries of the loop integration domain selected by the momentum mapping in Eq.(3.9) We can now express the phase space factorization as dΦ m q 1δ (q 1 ) R 1 (ξ 1 , v 1 ) = µ 2ε dΦ m+1 θ(y bc − y ac )J(y ac , y bc )ξ 1 (y ac , y bc ) . where dΦ m (dΦ m+1 ) is the customary differential phase space element for m (m+1) partons in the final state.
In the region of the real phase space where s bc < s ac (i.e. the one where p b and p c can be collinear), the dual contributions in which q 2 (green propagator in Fig.(3a)) is set on-shell can be used to match the soft and the collinear singularities [46]. To this end, a second mapping can be used with all the other partons left unaltered. In terms of (y ac , y bc ) and (ξ 2 , v 2 ), we have and the associated jacobian is (3.20) The inverse mapping used in [46] leads to which gives us the definition domain through Finally, phase space factorization is given by Given a pair (i, j) of partons with momenta p i and p j , LTD reconstructs the entire triangle in Fig.(3a) by cutting the propagators of q 1 , q 2 and q 3 . Nonetheless, to match single soft and collinear singularities associated with the (a, c, b) real parton phase space configurations, we just need to collect the dual contributions with q 1 and q 2 set on-shell, selecting the loop integration domains by using the functions R 1 and R 2 defined in Eqs. (3.16) and (3.23) [46]. In the construction of a working subtraction scheme, we may follow a slightly different approach, though. In fact, the region where s bc < s ac can be covered by simply exchanging p a ↔ p b and p i ↔ p j in the momentum mapping of Eq. (3.9). Then, to cancel the corresponding singularities, we just need to exchange p i ↔ p j in the dual counterterms with q 1 on-shell. This corresponds to cutting the propagator of q 1 in the triangle of Fig.(3b) instead of the propagator of q 2 in the triangle of Fig.(3a). In this way, we always put q 1 onshell, avoiding to introduce different counterterms, and we use one mapping instead of two. Moreover, if we chose to set on-shell q 2 and use the momentum mapping in Eq. (3.18), this would reflect into a more involved counting of the counterterm matching the singularities of the (m + 1)-parton phase space. More details on this subject are given in Appendix A.

Singular behavior of real amplitudes
In this Section we analyze the factorization properties of tree-level amplitudes with m + 1 final state partons in the soft and collinear regions of the phase space [5,63,64]. The formulas of this section will be our reference point to check that the dual counterterms we are going to define have the same local behavior of the real cross section in the IR singular regions.
The singular behaviour of an (m + 1)-parton tree-level amplitude does not depend on the details of its structure. Furthermore, it turns out that in the soft [64] and collinear [63] regions the (m + 1)-parton amplitude is factorizable with respect to the m-parton one.
Let us first consider the limit where two momenta, say p a and p c , become collinear. By definition, this configuration is parametrized by where p and n are light-like vector (with p individuating the collinear direction), k ⊥ (k 2 ⊥ < 0) is the transverse component (k ⊥ · p = k ⊥ · n = 0) and z is the fraction of momentum p carried by p a in the collinear limit which is approached when k ⊥ → 0.
In the collinear limit, the (m + 1) tree-level matrix element squared has the following behaviour [63] where the m-parton matrix element on the right-hand side is obtained by replacing the partons a and c in the (m + 1)-parton matrix element by a single parton denoted by ac (the emitter), with momentum equal to p a + p c and other quantum numbers (flavour, colour) given by the following rule: anything + gluon gives anything, and quark + antiquark gives gluon. The matrix in color and spinP (ac),c on the right-hand side of Eq.(3.26) is the treelevel d-dimensional Altarelli-Parisi splitting function. It is a matrix acting on the spin indices of the parton ac in m 1, . . . , m + 1| and |1, . . . , m + 1 m . Its value depends on whether the partons ac and a are quarks, antiquarks or gluons, according to where s| · |s and µ| · |ν denote the components of the splitting functionP (ac),a in the spin space of the parton ac. In Eqs. (3.27), the label g stands for gluon, while q stands both for quark and antiquark.
In case the parton a belongs to the initial state, we use a quite different parametrization describing the splitting process a → ac + c, where c is a final state parton and ac is the parton entering the hard scattering, whose quantum numbers are assigned according to the following rule: if a and c are partons of the same type, then ac is a gluon, while if a is a fermion (gluon) and c is a gluon (fermion), then ac is a fermion (antifermion). For the case of emission from an initial state parton, in the limit k ⊥ → 0, the tree-level matrix element behaves as follows where the m-parton matrix element on the right-hand side is obtained by replacing the partons a and c in the (m + 1)-particle matrix element by the initial state parton ac, which takes the role of a. The tree-level Altarelli-Parisi splitting functions for the case of initial state radiation can be obtained from the ones in Eqs. with F (q) = F (q) = 1, F (g) = 0. They are given by The limit where a gluon, say p c , becomes soft, is parametrized by where p µ is an arbitrary vector that specifies the direction along which the parton c approaches the soft limit. When λ → 0, we have the following behaviour for the (m+1)-parton matrix element squared [64] where the m-parton matrix element on the right-hand side is obtained by removing the soft gluon p c from the (m + 1)-parton matrix element. As we can see, the m-parton matrix element is not exactly factorized, because of colour correlations induced by the colourcharge operator T a T b . Nonetheless, we still have exact factorization in processes with m = 2 and m = 3 since, in these cases, the product of two colour-charge operators is always proportional to a linear combination of Casimir operators. If there are partons in the initial state, Eq.(3.33) still holds: we just need to extend the sum over more pairs, also including the initial state partons.
Let us now move to the case of massive partons in the final state [65,66]. If two momenta become collinear, and at least one of them is massive, there is no collinear singularity. However, the matrix element squared is enhanced for very small values of the parton masses. To take this effect into account, we parametrize the momenta p a and p c of the two collinear partons in the following way [6] where m a , m c and m ac are the masses of the partons a, c and the emitter parton ac, respectively, p is a time-like vector which individuates the collinear direction, while n, k ⊥ and z have the same meaning as in the massless case (see Eqs. (3.25)). We then consider the behavior of the tree-level matrix element under the following uniform rescaling k ⊥ → λk ⊥ , m a → λm a , m c → λm c m ac → λm ac (3.35) in the limit λ → 0. We have [6] m+1 1, . . . , m + 1|1, . . . , m where the m-parton matrix element on the right-hand side is obtained in the exact same way as for Eq. (3.26). The functionP (ac),a on the right-hand side of Eq.(3.36) is the generalization of the d-dimensional Altarelli-Parisi splitting function to the massive case (the dependence on the masses m a , m c and m ac being denoted by {m}). As in the massless case, it is a matrix acting on the spin indices of the parton ac in m 1, . . . , m + 1| and |1, . . . , m + 1 m , whose expression depends on the type of partons involved in the splitting process ac → a+c. Denoting quarks and antiquarks by q and gluons by g, we have where, as usual, s| · |s and µ| · |ν label the components of the splitting functionP (ac),a in the spin space of the emitter parton ac. The functionP g,g related to the splitting process g → g + g remains the same of Eq.(3.27d), since gluons are massless in any case. Unlike the collinear limit, the soft one leads to a real singularity in the tree-level matrix element squared, whatever value the emitter mass has. The soft emission can be parametrized as in the massless case (see Eq.(3.32)). The (m + 1)-parton matrix element squared then behaves as follows where the m-parton matrix element on the right-hand side is obtained, as in the massless case, by removing the soft gluon p c from the (m + 1)-parton matrix element. The difference between Eq.(3.38) and its massless analogue (Eq.(3.33)) lies in the presence of terms proportional to the mass of the particles that radiate the soft gluon.

Dual Subtractions
Selecting a pair of hard partons (i, j), we can apply LTD to extract the associated IR divergences by cutting the q 1 and q 2 propagators in Fig.(3a). As anticipated in section 3.1, in place of the dual contributions with q 2 set on-shell, we can exchange p i ↔ p j (see Fig.(3b)) and describe the rest of the divergences again by cutting q 1 . Nonetheless, for completeness in the following we will present dual subtraction formulas including also the cut over q 2 .
Contributions of the cut over q 1 Let us start by considering the case in which a gluon with momentum p c is radiated collinear to a quark of momentum p a and absorbed by a spectator of momentum p b . The corresponding dual subtraction term is constructed by summing the dual representation of two contributions: • the triangle involving the emitter and the spectator in the final state, as well as the bubble depending on the emitter momentum, that we collectively call V (1) ac,b ; • the wave-function renormalization of the emitter, that we denote by G Let p i , p j and q 1 be the virtual sector momenta of the emitter, the spectator and the virtual radiation, respectively, which are matched to the real sector momenta by the mapping in Eq.(3.9). The dual counterterm is then given by where N in includes all the non-QCD factors, S {m} is the Bose symmetry factor for identical partons in the final state, |1, . . . , m m is the m-parton Born amplitude and T ac T b is a colour-charge operator acting on the colour indices of the partons i and j, respectively. We have used the labels ac and b for these operators because the flavours of the partons i and j are the same of the partons ac and b of the real sector, respectively. The functions V qg,b groups together the dual contributions from the reduction of the virtual amplitude. In particular, in V (1) qg,b we can identify the dual representation of the following scalar and we have extracted only terms where the internal momentum q 1 is set on-shell in the general formula for LTD in Eq.(3.5). Moreover, we have not included the B 0 (p j ) bubble integral neither the wave-function renormalization of the spectator parton because none of them develop a divergence cutting the q 1 internal propagator and integrating over the phase-space region selected by R 1 . As stated above, G qg,b comes from the quark wave-function renormalization [46,48] (3.42) with M being the on-shell fermion mass that we assume to be zero for the time being, postponing to Section 4 the case of radiation off massive fermions. Note that the integrand of Eq.(3.42) depends on an auxiliary momentum p, but the integral does not. When contracting ∆Z quark with the Born matrix elements, a dependence on the spectator colour-charge operator T spec ≡ T b can be introduced by expressing the Casimir C F = −T 2 emit = −T 2 ac according to the colour conservation relation Then, in each term of the sum over the different spectators we can choose the auxiliary momentum p to coincide with the spectator momentum p j . In this way, the wave-function renormalization counterterm is decomposed as By applying LTD, selecting the dual contribution with q 1 on-shell and restricting the integration domain to the R 1 region, the definition for the contribution G (1) in Eqs. (3.39) and (3.40) follows. Note that, the presence of the spectator momentum and the restriction of the integration domain to R 1 , introduces a dependence on the s ij invariant. However, such a dependence does not affect the singular behaviour and only reflects in the finite part of the integral.
In case a gluon is radiated collinear to another gluon, the corresponding dual counterterm σ where the symmetrization p a ↔ p c has to be performed in the right-hand side of the momentum mapping in Eq.(3.9) and gives rise to a distinct counter event. Once again, in Eq.(3.45) we have extracted only the dual contributions with q 1 on-shell, obtaining In Eq. (1)µν gg,q , impedes the factorization of the Born amplitude squared. Nonetheless, this is a good feature for the dual counterterm, since Eq.(3.27d) shows us that the same kind of Lorentz structure appears in the (m + 1)-parton matrix element squared when we approach the limit of a gluon being emitted collinear to another gluon. Note however that, in the Feynman gauge, all the terms in Eqs. (3.46) proportional to −g µν can be easily contracted by results from the same scalar integrals of Eqs.(3.41), once we have embedded −g µν into the Born matrix elements as in Eq. (3.47). Note that the coefficient of the bubble contribution in Eqs. (3.46) is decreased by a factor 2 with respect to the case of a gluon emitted from a quark, as a result of the virtual amplitude reduction. The counterterm G (1)µν gg,b , on the other hand, is obtained by application of LTD to an integrand representation of the gluon and ghost contributions to the gluon wave-function renormalization counterterm G µν (3.48) The above expression is also justified by both its ε-poles structure and its tensorial properties. Given that the integrand representation G µν gg,b of the gluon wave-function renormalization is scaleless, its unconstrained integration over the loop momentum vanishes.

The integration of G
(1)µν gg,b (performed in Section 3.5, Eq.(3.102)) exhibits a single IR pole that coincides, up to a minus sign, with the correct UV pole for the gluon wave-function renormalization. As for the tensorial structure of the counterterm, note that, once defined we have N (1)µν p 1µ = N (1)µν p 1ν = 0, so that the contraction between the integrand on the right-hand side of Eq.(3.46b) and the tensor p 1µ p 1ν is equal to zero. Moreover, it is interesting to compare the tensor N (1)µν with the one used in [5] for the dipole associated with the emission of a gluon from another gluon, that we report here for convenience In fact, by using the mapping in Eq.(3.10) to express the vector w µ in terms of the virtual sector variables, we have Since, by gauge invariance, we have Eq.(3.51) tells us that the two tensors N (1)µν and M µν differ only by gauge terms. The dual subtraction for the remaining case of a collinear quark-antiquark pair includes the quark contribution to the gluon wave-function renormalization, and takes no contribution from the loop correction. We have which provides an integrand representation for the renormalization counterterm mentioned above. In a way analogous to the gluon and ghost contribution, both the pole and the tensorial structure justify its expression. Note that the insertion of the colour-charge operator T ac T b has been forced, since the sum of Eq.(3.43) simplify the Casimir C A at the denominator of Eq.(5.16). However, as it happens for the other renormalization counterterms, the dependence on the spectator plays no role in the singular behaviour. For completeness, we show the integrand representation of the whole gluon wavefunction renormalization counterterm, obtained by combining gluon, ghost and quark contributions with p an auxiliary momentum. Finally, note that the roles of quarks and antiquarks in the dual counterterms defined in this Section are interchangeable: when a quark takes the place of an antiquark or vice versa, we just need to change the sign of their momenta in the wave-function renormalization counterterms.

Contributions of the cut over q 2
The dual counterterms obtained by selecting the terms where q 2 is set on-shell by LTD are for the emission of a gluon from an antiquark as well as for the emission of a gluon from another gluon and (3.62) for a collinear quark-antiquark pair. Note that, with the notation σ D (k) ac,b , we indicate that q k is set on-shell, c is the radiated parton and parton b is the spectator if k = 1, while this role is played by a if k = 2. Analogous considerations hold also for V

Singular behaviour of the dual counterterms
In this Section we prove that the dual counterterms defined in Section 3.3 locally cancel the corresponding singularities of the (m + 1)-parton tree-level cross section. We do this by applying the momentum mappings of Section 3.1 and analyzing the singular behaviour of the dual counterterms in the infrared and collinear limits of the (m + 1)-parton phase space. We will see that, in the singular regions, the local behaviour matches exactly the one of the real matrix element squared shown in Section 3.2.
In the following, we denote by q 1 -cut method the strategy in which, to collect the singularities associated to the emission from the pair of hard partons i and j, we cut only the propagator connecting them (q 1 in Fig.(3)). Furthermore, we denote by q 1 -q 2 -cut method the method in which both the cut contributions over q 1 and q 2 are included to collect the above mentioned singularities. Nevertheless, we anticipate here that the q 1 -cut method will turn out to be more convenient.
Let us start by analyzing the region where the momentum of a gluon, say p c , becomes soft. If we substitute the parametrization of Eq.(3.32) into the right-hand side of Eq.(3.9), we obtain the following behaviour in the limit λ → 0 Since in what follows we are going to use the phase space factorization of Eq.(3.17), it is useful to evaluate the behaviour of the product J 1 ξ 1 . This can be done by substituting the parametrization of Eq.(3.32) into the right-hand side of Eqs.(3.13) and (3.14), obtaining integrals over Φ m+1 , we get the following singular behaviour for the corresponding dual cross sections . Note that we have used T ac = T a because c is a gluon. In order to obtain the total dual cross section for the soft emission of the gluon c, we need to perform a sum over all the possible emitter-spectator pairs and multiply times a factor 1/(m g + 1) that turns S {m} into S {m+1} (see Appendix A). The result is the following Note that, in the limit λ → 0, the matrix element |1, . . . , m m can be considered as obtained from the (m + 1)-parton one by removing p c , since in this limit we have Eq.(3.66) exhibits the same soft behaviour of the (m + 1)-parton matrix element squared in Eq. (3.33). Therefore, the local cancellation of soft singularities between the dual and the real cross section has been proved.
Let us now move to the collinear limit. Consider the case of a parton with momentum p c emitted collinear to a parton of momentum p a and absorbed by a spectator of momentum p b . The collinear limit is parametrized in Eqs. (3.25) which, once plugged into the right-hand side of Eq.(3.9), leads us to In Eq.(3.67b), we have used the fact that, in the matrix elements, p i → p and, consequently, m 1, . . . , m| p µ → 0 as well as p ν |1, . . . , m m → 0 because of gauge invariance. Next, we need the behaviour of the product J 1 ξ 1 in the collinear limit. This can be obtained by inserting the collinear limit parametrization (3.25) into the right-hand sides of Eqs. (3.13) and (3.14). We have The behaviour of the total dual cross section for a given splitting process is obtained by performing a sum over all the possible spectators, multiplying σ gg,b times a factor 1/(m g + 1) and σ Since the only dependence of the integrands of Eqs.(3.69) and (3.70) on the spectator lies in the colour-charge operator T spec ≡ T b , we can easily perform the sum over the spectators by using the colour conservation relation spec, spec = emit In this way, using T 2 q = C F and T 2 g = C A , we get By comparison with Eq.(3.26), we see that the integrands of Eqs.(3.73) share the same local behaviour of the (m + 1)-parton amplitude squared. Therefore, we have proved also the local cancellation of collinear singularities between the dual and the real cross section.
In case dual counterterms with q 2 on-shell are considered, to study their soft behaviour we need to substitute the parametrization of Eq.(3.32) into the mapping of Eq.(3.18) instead of Eq.(3.9). We have while, using Eqs. (3.19) and (3.20), we get aq,q ∼ O(λ −1 ). We have used T bc = T b because c is a gluon. The combination of soft contributions given by cutting over q 1 , Eq.(3.65), and q 2 , Eqs.(3.76) reads More details on the matching among the singularities of the real matrix elements with the dual counterterms are provided in Appendix A.
If parton b and c become collinear, the singular behaviour of the dual counterterms with q 2 on-shell can be obtained by replacing p a → p b in the parametrization of Eqs. (3.25) and substituting the result into the momentum mapping of Eq.(3.18). We have (3.79b) Since p j → p, we have used again that m 1, . . . , m| p µ → 0 and p ν |1, . . . , m m → 0 for gauge invariance. To obtain the collinear behaviour of the product J 2 ξ 2 , we substitute the parametrization of Eqs.(3.25) into Eqs. (3.19) and (3.20). We have Proceeding in a way analogous to the case of counterterms with q 1 on-shell, we obtain The sum of the dual counterterms with q 1 and q 2 on-shell over all possible spectators, gives ac,b are defined in Eq.(3.78)) so reproducing the correct p b p c collinear limit.

Integrated Dual Subtractions
In this section we show how to integrate the dual counterterm over the (constrained) loop momentum. The result lives on the m-particle phase space and, therefore, can be integrated together with the virtual cross section.
Let us start with the case of a quark or an antiquark as the emitter. In terms of the dimensionless variables (ξ 1 , v 1 ), the subtraction terms in Eq.(3.40) are given by where the measure [dξ k dv k ] is defined as The integrands in Eqs. (3.85) have been expressed in terms of the loop momentum variables (ξ 1 , v 1 ), defined in the center-of-mass frame of p i + p j by Eq.(3.12), using .87) and (removing the subscript k for easy of notation) where, in the square bracket of Eq.(3.89a), we have kept in separated round brackets the contributions of the scalar triangle (former) and bubble (latter) integrals. By substitution into Eq.(3.39), we get which has the same poles of the terms in Eq.(2.6) with γ i = γ q and T i = T q , thus confirming that we have extracted in the proper way the IR singular behaviour from the virtual cross section. gg,b of Eq.(3.40) times −g µν . The only difference lies in the relative coefficient of the scalar bubble with respect to the scalar triangle, which is multiplied by a factor 1/2 in the case of a gluon as the emitter, as we have already discussed. Therefore, using Eqs.(3.86), (3.87) and (3.88), we have and the result of the integral is where the remaining tensor −g µν can be contracted with the spin polarization indices of the parton p j in the m-particle matrix elements, leading to m 1, . . . , m|T i T j |1, . . . , m m as in Eq. (3.47).
The integration of the non-trivial tensor structure in G (1)µν gg,b requires particular attention. Let us focus on the term proportional to q µ 1 q ν 1 in Eq.(3.46b) By Lorentz covariance, the result of the integral must be of the type since −T µν A 00 g µν = A 00 and T µν p µ i = T µν p ν i = T µν p µ j = T µν p ν j = 0. The coefficient A 22 can be obtained by contraction with the tensor p iµ p iν /(p i · p j ) 2 . Therefore we have Now consider the integrals which represent the terms in Eq.(3.46b) proportional to q µ 1 p ν j and p µ j q ν 1 , respectively. By Lorentz covariance, these integrals must be of the type Again, because of Eq.(3.52), the only non-vanishing contributions are B 2 p µ j p ν j and C 2 p µ j p ν j . These can be extracted by contraction with p µ i /p i · p j and p ν i /p i · p j , respectively. Therefore, we can write Note, by looking at Eqs.(3.96) and (3.99), that the term proportional to p µ j p ν j in I µν 00 + I µν 01 + I µν 10 is opposite to the one proportional to p µ j p ν j in Eq.(3.46b). Therefore, we can write which, in terms of the dimensionless variables (ξ 1 , v 1 ), becomes The result of the integral is the following where, again, the tensor −g µν can be contracted with the m-particle matrix elements, leading to m 1, . . . , m|T i T j |1, . . . , m m . We now have to insert Eqs.(3.92) and (3.102) into Eq.(3.45), contracting −g µν as in Eq. (3.47). Note that we do not have to consider the symmetrization p a ↔ p c , since it is already taken into account in the virtual sector. We obtain Following the same reasoning, the integration of the quark contribution to the gluon wave-function renormalization leads to The sum of the integrated dual cross sections in Eqs.(3.103) and (3.105) has the same poles of the terms in the sum of Eq.(2.6) with γ i = γ g and T i = T g , as we expected. The dual counterterms with q 2 on-shell can be integrated in a way analogous to the case of q 1 on-shell. The results are the following It is interesting to note that, if we use only terms with q 1 on-shell, we are able to entirely reconstruct the real part of the triangle associated with p i and p j in Eq.(2.4), with no remnant left out. In fact, as we can appreciate in Eqs.(3.89a) and (3.92), the dual contributions coming from the scalar triangle in the sum V (1) bg,a lead us to It should be noted that this result depends on the function R 1 used to select the loop integration domain which, in turn, is dictated by the mapping. Therefore, the reconstruction of the whole triangular function may be considered an accident. However, this could also suggests a connection between the particular mapping of Eq.(3.9) and the structure of the three-point scalar function. From now on, we will use only the q 1 -method.

Masses in the final state
The presence of massive quarks in the final state does not alter the reasoning of section 2. The difference is that now we have to apply the TLD theorem with massive propagators and that we need a generalization of the mapping in Eq.(3.9) to take into account the mass of the fermions. For illustrative purposes, in this section we limit ourselves to the case of a single emitter-spectator pair of fermion-antifermion with the same mass. This example highlights the main differences with the massless case and can be easily generalized to the cases of pairs with different masses or one massive and one massless particle.

Mapping between virtual and real sector
Let us denote by M the common mass of the quark and the antiquark. The mapping in Eq.(3.9) can be generalized as follows where α abc , β abc and γ abc are defined by As in the massless case, the mapping in Eq.(4.1) automatically verifies momentum conservation p a + p b + p c = p i + p j and on-shell conditions p i 2 = p j 2 = M 2 . We also choose again the region p a · p c < p b · p c to be the definition domain of the transformation in Eq.(4.1). Note that, in the massless limit (M → 0), α abc → 0, β abc → 1 and γ abc → 1 − 2(p a · p c )/s abc . Therefore, in this limit, the mapping in Eq.(4.1) reduces to the one in Eq.(3.9) that we have used for the massless case.
Working out the q 1 -method, we set on-shell only the loop momentum flowing into a massless propagator. For this, we can again use the parametrization of Eq.(3.12) to assign q 1 in terms of the dimensionless variables (ξ 1 , v 1 ) in the center-of-mass frame of p i + p j .
In order to obtain the phase space factorization, we may want to express the momentum mapping in Eq.(4.1) as a relation between the kinematic invariants (y ac , y bc ) and the variables (ξ 1 , v 1 ). The momentum mapping in Eq.(4.1) implies ξ 1 = y ac + y bc v 1 = y ac − (y ac + y bc )(y ac + α abc ) (y ac + y bc )(1 − y ac − 2α abc ) where α abc can be easily expressed as a function of (y ac , y bc ) directly from its definition in Eq.(4.2).
To compute its expression in terms of (ξ 1 , v 1 ), we start by considering the inverse of the momentum mapping in Eq.(4.1), that is given by [48] p c = q 1 wherep i andp j are massless momenta related to p i and p j by The expressions of α ij = α abc , β ij = β abc and γ ij = γ abc in terms of the virtual sector variables are the following

Inverting (4.3) we get
that can be used, together with (4.8), to express R 1 in terms of (ξ 1 , v 1 ), obtaining (4.10) In the region y bc < y ac the prescription of the q 1 -cut-method consists into exchanging both p a ↔ p b and p i ↔ p j , i.e. we consider the diagram in Fig.(3b) and we set again q 1 on-shell.

Dual counterterms and singular behaviour
In the case of a quark-antiquark massive pair as emitter and spectator, the application of LTD to the triangle and the bubble contributions from the virtual cross section leads to the same result as for the massless case. In fact, according to Eq.(3.6), the relevant dual propagators are given by which do not depend on the mass of the fermions. The same does not happens if we consider the dual contributions with q 2 on-shell (q 1 -q 2 -method ).
With this in mind, we define the dual subtraction term as Here, the first contribution to V (1) qg,q represents the triangle scalar function coming from the reduction of the virtual amplitude. Contrary to the massless case, the coefficient of the bubble integral of the external massive momenta coming from the reduction, is highly nontrivial. Nevertheless, this scalar integral is IR finite. However, when the mass of the fermion is very small relative to the other scales of the process, this contribution exhibit a logarithmic enhancement in both the real and the virtual sector. It is then convenient to keep the bubble in the counterterm V (1) qg,q with the coefficient of the massless case. The dual counterterms G 1,qg,q and G (1) 2,qg,q come from the integrand representation of the wave-function renormalization of the emitter and the spectator, respectively, reported in Eq.(3.42) [48].
The procedure to prove that the dual subtraction in Eq.(4.12) has the same singular behaviour of the real amplitude is algebraically challenging but straightforward. For this reason, here we limit ourselves to just list the steps of the proof. As for the quasi-collinear limit, the first thing to do is to express the scalar products appearing in Eq.(4.13) in terms of p a · p b , p a · p c , p b · p c and M 2 by using the momentum mapping in Eq.(4.1). Then, in order to test the collinear behaviour, the invariants p a · p b and p a · p c have to be written as a function of z, k ⊥ , p, M 2 and n via the parametrization in Eqs. (3.34). At this point, the uniform rescaling in Eq.(3.35) can be performed and the dual subtraction in Eq.(4.12) can be expanded in series of λ. What is found is that the leading order perfectly matches the right-hand side of Eq.(3.36) withP (q,q) (z, k ⊥ , {m}, ε). In the soft region we need to use Eq.(3.32) instead of Eqs. (3.34), and then to consider again the series of λ. The result obtained with this procedure is in agreement with the right-hand side of Eq.(3.38).

Integrated dual counterterms
As in the massless case, we need to integrate the dual counterterms over the loop momentum, in order to obtain a result which lives on the m-particle phase space and is ready to cancel the poles of the virtual amplitude.
First of all, the massive momenta p i and p j are given, in their center-of-mass frame, by the following expressions (4.14) The scalar products between the loop momentum q 1 parametrized in Eq.(3.12) and the external momenta of Eq.(4.14) turn out to be These integrals can be computed analytically. Their poles are given by By inserting Eq.(4.16) into Eq.(4.12) we obtain, for the dual cross section, the following pole which is doubled when we add σ qg,q . The pole structure just obtained matches the contribution to the one-loop amplitude associated with a massive quark-antiquark emitterspectator pair [65].

Initial state radiation
It is well known [63] that in presence of partons in the initial state the singularities of the virtual matrix element do not cancel the ones of the real contribution and the proper definition of the NLO cross section requires the subtraction of the remaining collinear initial state singularities. The reasoning of Section 2 which allows to reconstruct the universal IR behaviour of the renormalized one-loop amplitude holds also for the case of partons in the initial state. So that, we will start by the singularities of the virtual matrix element knowing already that they will cancel all but the collinear initial state singularities of the real matrix element. We start by considering a quark emitter and a gluon radiation. Following the momentum assignment of the left panel in Fig.(4) (note that now momenta p i and p j are incoming so that in the expressions below there are different signs with respect to the expressions reported in Section 3) and recalling the same steps of Section 3 the singularities of virtual contribution associated to the emission of a virtual gluon by an hard quark carrying momentum p i absorbed by the hard parton carrying momentum p j is given and R 1 is any function that cuts UV divergences, selects the region q 1 p i and leaves out the one where q 1 p j . Let's focus on the case in which partons i and j are in the initial state. We will extend the results to the other cases in a moment. When, in the previous sections, we have studied the case of final state radiation, we enforced the constraint of preserving the invariant mass of the Born system. This time at Born level we are producing a final state with invariant mass s ij and so enforcing the conservation of the invariant mass when we cut the scalar loop functions would generate by definition a phase space configuration sitting on the soft limit for the emitted gluon q 1 . We now consider the real phase space volume corresponding to an initial state configuration with two partons with momenta p a and p b such that s ab > s ij . In this case one can exploit the phase space factorization to map the real phase space onto the Born phase space (with initial momenta p i and p j ) and, following the same approach of the previous sections, identify the (extended) loop momentum with the momentum of the radiation. It is a straightforward exercise to show that if parton a is a quark, the collinear divergences of the real matrix element associated to a gluon emission can be described in terms of the above formula multiplied by the momentum fraction x of the relation p i = x p a . We anticipate that this extra factor of x has to be included for all the cases with an initial emitter (i.e. all initial-initial and initial-final cases). To show that the dual counterterm (including this factor of x) cancels the initial state collinear singularity of the real matrix element, we adopt the momentum mapping for initial-initial singularities of ref. [5]. Apart from the boost on the Born system, the mapped momenta are given by where p c is the momentum of the radiated gluon in the real phase space configuration. We have then and expressing all the scalar products in terms of the invariant built with the momenta of the collinear pair (p a · p c ) and the variables x and v, we have Substituting the scalar products in the right hand side of Eq.(5.1) and multiplying by x we find which in the limit v → 0 coincides with the first of Eqs. (3.31). This counterterm has to be integrated at fixed s ab and imposing the condition p a · p c < p b · p c we constraint the integration region to As usual, the integration has to be performed in the sense of distributions and gives rise to a distribution in x. Finally, to compensate for the different flux factor of the real contribution (s ab ) with respect to the one of the virtual contribution s ij = x s ab , for the analytic integration we multiply the counterterm times x. Also the cases of initial-final and final-initial singularities will be addressed with the corresponding momentum mapping proposed in ref. [5] that enforces one rescaled momentum in the initial state. For this reason also for these cases the integrated version of the counterterms is multiplied times this compensating factor of x. A small difference with respect to the case of initial-initial is the presence of an extra jacobian factor that is absent (i.e. 1) for initial-initial, but is different from 1 for the other two cases, initial-final and final-initial. The result for the integrated version of the dual counterterm in Eqs. (3.31) and for all the others of the massless case discussed below are collected in Appendix B.
Let us now consider the case of a gluon emitting a gluon in the initial state. When we discussed the case of final state emission in Section 3.3 we introduced the dual counterterm that we report here with incoming momenta for convenience After we introduced the mapping that splits the emitter and spectator momenta p i and p j into the collinear pair p a , p c and the new spectator p b we proved that the collinear limit of the two partons p a and p c in real matrix element is reproduced by this virtual contribution by summing over the two possible final state gluon configurations obtained by exchanging the momenta p a and p c of the collinear pair. Such symmetric role of the two radiated partons is lost in the case of radiation from the initial state (note that in the equation above we are not doing the symmetrization), but still of course the above formula correctly reproduces the poles of the virtual matrix element. Nevertheless, if we multiply this expression times x and consider the extended kinematic range for the loop momentum associated to the possible radiation in a real phase space kinematic with s ab > s ij this time we cannot neglect the contribution from the other cut gluon propagator in the right panel of Fig.(4). Of course, all three kinds of integrals, not just the scalar triangle in the figure, but also the scalar bubble and the wave-function renormalization constant have a second cut contribution, although their sum has to vanish in the soft limit x → 1. The presence of a square denominator in the wave-function renormalization is not particularly problematic per se thanks to the general results of ref. [46] on the computation of dual integrals with denominators with powers greater than 1. However, to simplify the construction here we include these three extra contributions expressing them in terms of the same integration variables of the others. To this aim we invert the loop momentum and shift it so that the cut propagator always has momentum q 1 (see right panel of Fig.(4)). In this way, the second cut contributions from the scalar bubble and from the wave-function renormalization constant turns out to be identical to the ones brought by their first cut (as one could naively expect) while a genuine new contribution comes from the scalar triangle. The dual subtraction is then given by To study the collinear limit of the expression above we start by analyzing the non trivial tensor part. The vector is orthogonal to the emitter momentum p i so that in general it is a superposition of k µ ⊥ (that is not massless) and p i (that is massless and orthogonal to k ⊥ ) given by The tensorial components proportional to p µ i and p ν i are gauge terms and will not contribute so that multiplying and dividing the tensor part by N 2 the counterterm in the square breacket of Eq.(5.9) can be written as Inserting the scalar products in Eqs.(5.4) and (5.5) and taking the v → 0 limit, the expression above times x coincides with the collinear limit in Eqs.(3.31) for a gluon emitting a gluon. We now move to the counterterms for the splittings in which the initial parton is different from the hard (emitter) parton. We first discuss the case of a quark emitting a quark transforming into an hard gluon for the underlying Born configuration. In this case the dual counterterm is built from the quark contribution to the gluon wave-function renormalization constant and is given by (including the factor of x for the case of an initial emitter) (5.16) Following the same steps as for the case of a gluon emitting a gluon it is easy to show that in the limit v → 0 this counterterm reproduces the functional form of the third equation of Eqs. (3.31). The last case is when an initial gluon radiates an antiquark transforming into a hard quark. In this case we find that the sole contribution of the quark wave-function renormalization constant (times x) is not sufficient to get the correct collinear limit of the real matrix element. We have the correct limit including a contribution given by minus the second cut of the virtual part associated to a gluon emitter. The result is: Once the scalar products are substituted in the formulae above, in the limit v → 0 the resulting expression reproduces proper collinear behaviour of Eqs. (3.31). We stress that the dual dipoles for the final-final case describe also the final-initial case after the corresponding jacobian factor is included to translate the radiation variables into the loop momentum.
In the same way the formulae of this section for the initial-initial case describe also the initial-final case once the proper jacobian factor is used. We have then completed the list of all possible counterterms. A collection of all the relevant formulas for all the countertems and their integrated version is given in Appendix B.

Applications
In this Section we will show a small selection of examples where we apply the construction presented in the previous sections. For illustrative purposes, in the present section we limit the discussion to a set of relatively simple processes. In each case, a comparison is made among the results obtained using dual counterterms and the ones obtained using Catani-Seymour dipoles.

γ * → 2 jets at NLO
At the lowest order, the two-jet production in e + e − collisions has just a quark-antiquark pair in the final state with momenta p 1 and p 2 , respectively. In the real sector, we have a single sub-process where an additional gluon is radiated with momentum p 3 , while the quark and the antiquark have momenta p 1 and p 2 , respectively. Therefore, the only emitterspectator pair is the one composed by the quark and the antiquark, with the emitter role being played once by the quark and once by the antiquark. This means that we need two dual counterterms: σ qg,q for the case where the emitter is the quark and σ qg,q when that role is played by the antiquark.
In the Monte Carlo implementation, once a real phase space configuration (p 1 , p 2 , p 3 ) is generated, the kinematic invariants s 13 and s 23 are compared. If s 13 < s 23 (R 1 region), the integrand of the dual cross section σ qg,q is activated. In any case then, the method utilizes a single counter event for each real phase space configuration.
In the virtual sector, we use the integrated dual counterterms of Section 3.5. In particular, for both counterterms the result of the integration is equal to the right-hand side of Eq.(3.90). Since there are only a quark and an antiquark in the virtual sector, we have T 1 T 2 = −C F . Therefore, we can write qq | 2 is the leading-order matrix element squared and we have neglected terms of order O(ε).
We now show the result for a differential prediction. We have chosen √ s = 125 GeV as the energy in the center-of-mass frame. The corresponding value of the strong coupling constant, α S = 0.11173799, has been obtained starting from α S (m Z ) = 0.117, m Z = 91.1876 GeV being the mass of the Z boson. The jet algorithm used in the computation is the Durham algorithm [67] with y cut = 0.1. Using this setup, in Fig.(5) we plot the differential cross section with respect to the rapidity of the most energetic jet. The excellent agreement observed in Fig.(5) among the two results validates our procedure and is found by plotting any other differential variable.

γ * → 3 jets at NLO
The basic sub-process contributing to the three-jet production in e + e − annihilation into a virtual photon is γ * → q(p 1 )q(p 2 ) g(p 3 ), while in the real sector one has the two subprocesses γ * → q(p 1 )q(p 2 ) g(p 3 ) g(p 4 ) and γ * → q(p 1 )q(p 2 ) q (p 3 )q (p 4 ), where the two quark-antiquark pairs may have equal or different flavours. For simplicity we will consider only the gluonic channel, where an additional gluon is radiated with respect of the virtual sector. The contribution of the other sub-processes can be accounted for by similar considerations.
By inspecting the virtual sector, we find six emitter-spectator pairs for the process γ * → q(p 1 )q(p 2 ) g(p 3 ), that are {qq, qg,qg} and the three pairs where the emitter and the spectator switch. The case of four quark production can be treated along the same lines. Once a real kinematic configuration (p 1 , p 2 , p 3 , p 4 ) is generated, the six kinematic invariants s 12 , s 13 , s 14 , s 23 , s 24 and s 34 are analyzed and six dual counterterms are activated. To evaluate them, six virtual configurations (p 1 , p 2 , p 3 , q 1 ) are built by using six times the inverse mapping in Eq.(3.9), every time with the proper set of real and virtual momenta.
To each emitter-spectator pair we associate the following dual counterterms terms of our algorithm.
In the virtual sector, we use the integrated dual counterterms of Section 3.5. Since there are three partons, a quark, an antiquark and a gluon, the colour algebra factorizes and one has T 1 T 2 = C A /2 − C F and T 1 T 3 = T 2 T 3 = −C A /2. Therefore, the integrated dual counterterms are given by where |A (0) qqg | 2 is the Born amplitude for the process γ * → qqg and we have neglected terms of order O(ε).
In Fig.(6) we show the differential cross section with respect to the rapidity of the most energetic jet using again √ s = 125 GeV as the energy in the center-of-mass frame, α S ( √ s) = 0.11173799 (starting from α S (m Z ) = 0.117) and the Durham jet algorithm with y cut = 0.05. Also in this case, we observe a perfect agreement with the same computation performed using Catani-Seymour dipoles. We refrain to include other plots because the same level of agreement is observed for any other differential distribution.

H → bb at NLO
In the present Section, we will make of use the formulas for the massive case derived in Section 4 to study the Higgs boson decay into bottom quarks. The virtual sector process is H → b(p 1 )b(p 2 ), while in the real sector an addition gluon is emitted, H → b(p 1 )b(p 2 ) g(p 3 ). Therefore, the counting of the dual counterterms is the same as for the process γ * → 2 jets analyzed in Section 6.1.
Once a real phase space configuration is generated, the proper counterterm in Eq.(4.12) is selected comparing the s 13 and s 23 invariants. As for the virtual contribution, we have renormalized the mass of the wave-function of the bottom quark on-shell. The virtual singularities are canceled by the integrated dual counterterms of Section 4.3 with T 1 T 2 = −C F . The differential distribution of the rapidity of the most energetic jet is plotted in Fig.(7). We have used m H = 125.09 GeV for the Higgs boson mass and α S (m H ) = 0.11263619 corresponding to α S (m Z ) = 0.118. The on-shell mass of the bottom quark has been set to m b = 4.78 GeV. The jet algorithm used in the computation is the Durham algorithm with y cut = 0.01. A fairly good agreement with the computation performed using Catani-Seymour dipoles is observed.

Drell-Yan pair production plus 0 or 1 jet at NLO
In this Section we consider proton-proton collisions where a W + boson is produced. We start by considering the Born process with no further radiation. In the virtual sector the basic sub-process is q(p 1 )q(p 2 ) → W + (p 3 ), while in the real sector we have q(p 1 )q(p 2 ) → W + (p 3 ) g(p 4 ). Since there are two partons in the initial state, we apply the method shown in Section 5. For the sake of comparison and for simplicity, we limit here to the quarkantiquark initiated sub-processes.
In the Monte Carlo implementation, we generate a real phase space configuration (p 3 , p 4 ) and compare the kinematic invariants s 14 and s 24 . If s 14 < s 24 , we apply the mapping in Eq.(5.3), with (p a , p b , p c ) = (p 1 , p 2 , p 3 ), to obtain a virtual configuration (p 1 , p 2 , q 1 ). Then, we evaluate the integrand of the dual counterterm σ D (1) qg,q in Eq.(5.7) and we add it to the real cross section. If s 14 > s 24 , we use the mapping in Eq.(5.3) with (p a , p b , p c ) = (p 2 , p 1 , p 3 ), instead. In the virtual sector we use the integrated dual counterterms of Section 5, with T 1 T 2 = −C F .
In Fig.(8) we show the differential cross section with respect to the rapidity of the W + boson. We consider proton-proton collisions at 8 TeV in the center-of-mass frame and we have used the MSTW2008NLO set of parton distribution functions with α S (m Z ) = 0.12018. The renormalization and the factorization scale have been set equal to µ R = µ F = µ = 80 GeV and we have used m W = 80.398 GeV for the W + boson mass, G F = 1.16639·10 − 5 GeV −2 for the Fermi constant and α S (µ) = 0.12264887 for the strong coupling constant. Also in this case, we observe an excellent agreement among the computation performed using dual subtractions and the one obtained with Catani-Seymour dipoles.
The case in which there is radiation also at Born level, beyond the initial-initial counterterms, requires the computation also of initial-final and final inital counterterms. For each possible emitter spectator pair the two kinematic invariants are computed and com- pared. The counterterm to subtract correspond to the smaller of the two invariants. In Fig.(9) we show the lepton and leading jet transverse momentum distributions respectively. Jets are formed using the antikt algorithm with resolution paramenter R = 0.5. Only the qq channel is included.

Higgs boson production in gluon fusion plus 0 or 1 jet at NLO
We consider the production of an Higgs boson in proton-proton collision in the framework of the effective field theory where the heavy top quark degree of freedom that mediates the interaction among the Higgs boson and gluons has been integrated out. At the Born and virtual sector the basic sub-process is then g(p 1 ) g(p 2 ) → H(p 3 ), while in the real sector we have g(p 1 ) g(p 2 ) → H(p 3 ) g(p 4 ). For simplicity, we limit here to a theory with only gluons, so that these are the only two partonic subprocesses for a NLO calculation. Since the Born process has a colourless final state and presents two partons in the initial state, we apply the method shown in Section 5. The mapping and the rest of the procedure is precisely the one already discussed for the case of the W boson production with the only difference that for the present case we use the formulas for a gluon emitting a gluon. In Fig.( Finally, we consider the process in which one hard radiation is present already at Born level. In this case the basic subprocess is g(p 1 ) g(p 2 ) → H(p 3 )g(p 4 ) and the real subprocess contains one further radiated gluon. Using the formulas for initial-initial, initial-final and final-initial gluon radiation from a gluon we cancel locally all the singularities obtaining fully differential predictions. In Fig.(11) we show the Higgs boson rapidity (left) and transverse momentum distributions (right) respectively. Again we observe a perfect agreement with the same computation performed using Catani-Seymour dipoles.

Conclusion
Starting from the analysis of the divergences of one-loop amplitudes reduced to scalar integrals, in the present paper we have built a subtraction formalism exploiting the LTD theorem. Beyond the singularities of the virtual diagrams, we have also to take into account the ones brought by the wave-function renormalization. An integrand-level representation for these counterterms suitable for massless and massive fermions can be found in [46] and [48], respectively. In the present work we provide an integrand-level representation of the gluon wave-function renormalization constant. Furthermore, our construction is naturally extended to the case of initial state singularities. Starting from the dual representation of the divergent part of virtual contribution, we cancelled the singularities of the real contribution  by simply feeding the former with the mapped momenta of the latter and multiplying times the momentum fraction x. We now expand a bit more on this last point. Let us consider a process in which a colourless final state is produced and, to further simplify the discussion, let's take the production of a pair or leptons from qq annihilation into a virtual photon, so ignoring the contribution of the Z boson. Normally, when calculating higher-order contributions in a proton-proton collision one starts from a leading order process producing a final state with a given invariant mass s ij . This might seems a complication in our case, but for the sake of the discussion, let's consider that our experiment is capable to see pairs of leptons only in a well defined narrow window of the invariant mass of the pair at the time. Once a value for the invariant mass is chosen we compute the Born, the virtual and the real contribution to the scattering process. It is well known that requiring a fixed final state invariant mass for the lepton pair the singularities will not cancel. Living on different phase spaces, s ij for the virtual and s ab > s ij for the real, they cannot match. Thanks to the remarkable property of factorization, this has been solved through a redefinition of the initial state in the form of a renormalization of the parton distribution functions. Let us now take a slightly different point of view. We could define the corrections, both real and virtual, to be initiated by the same center of mass energy s ab . The real part is precisely as before, but the virtual part is not contributing at all if s ab is greater then s ij . By the way, the divergences cancel if we consider s ab = s ij . The only difference now is that, if we have s ab > s ij and cut the loop amplitude in such a way to produce a gluon that is not soft, then the invariant mass of the leptons in the Born system will be lower then s ab and could be tuned to be the value s ij of the Born process. By doing so (and also multiplying times x) we observe the cancellation of the singularities. Of course, this way of cutting the loops produces a result that does not correspond to the original loop integral living on s ab , but it represents a generalization of the loop integral living on s ij for x < 1. In the present work we have exploited this feature that integrating such cut virtual loops, originally living on s ab , suffices to cancel both the loop divergences at s ij and the collinear initial state singularities. This fact provides a general strategy to use Loop-Tree Duality to compute cross sections with hadrons in the initial state. In fact our procedure can be applied also to the case of full amplitudes that are not reduced to scalar integrals. In practice, for the cases studied in the previous section, such as Drell-Yan and Higgs production in gluon fusion, the divergent part of the virtual contribution that we have extracted, differs very little (by finite terms and just one IR finite bubble) from the whole renormalized virtual contribution.
Acting such generalized virtual as a counterterm that cancels the IR divergences of the real part, there seems to be no room for an explicit presence of the factorization scale. Of course, we want to take benefit of parton evolution so the question is what is the factorization scale to be used in PDF's if we use a dual version of the (extended and cut) virtual contribution in the sense discussed above. This issue will be discussed elsewhere but we anticipate here a simple analytic strategy to answer to this question. It consists in deducing the corresponding factorization scale through the comparison among results obtained by performing the NLO computation in the two ways, using the cut extended loops on one side and the dipole subtractions plus the universal collinear counterterms (which bring the explicit presence of the factorization scale in the hard cross section) on the other.
On a technical ground, the structure of the dual counterterms is very close to the Catani-Seymour dipoles. In fact, even if we count a lower number of counter events per real phase space point with respect to the Catani-Seymour construction, in the latter it is possible to reduce the number of active dipoles by restricting the phase space region where they are computed to the neighbourhood of the singularities. On a more formal ground, the dual subtraction builds a direct link for the cancellation among virtual and real singularities.
We conclude by adding a few comments. The extension of the formalism presented here to the case of NNLO corrections is an interesting path that will be addressed elsewhere. We note that many ingredients are already available [68][69][70][71], including examples of full NNLO computation using LTD, momentum mappings and master integrals. Further applications of the work of the present paper would be the implementation of novel shower and matching schemes based on the dual subtractions discussed here. These subjects are beyond the aim of the present work and are left for future investigations.

A Counting of the dual counterterms
The computational scheme proposed in the current work exploits the dual representation of the loop amplitude with m partons in the final state in order to cancel soft and collinear singularities of the real matrix element living on the (m + 1)-parton phase space. The subtractions are extracted from the virtual sector on the base of the number of emitterspectator pairs in its m-parton phase space. Our starting point is represented by a set of dual counterterms C i,j (with i, j = 1, . . . , m and i = j), each associated with two functions V i,j and G i,j that depend on the emitter-spectator pair (p i , p j ). These are constructed by application of LTD to the scalar three-and two-point functions, C 0 (p i , p j ) and B 0 (p i ), respectively, and the wave-function renormalization counterterm ∆Z(p i ). In particular, C ij has the general structure in terms of the dual V and G The dual cross section is obtained by summing over all the possible emitter-spectator pairs where N in includes all the non-QCD factors, {m} indicates a sum over all the m-parton configurations and S {m} is the Bose symmetry factor for identical partons. The sum over the different emitters on the right-hand side of Eq.(A.2) can be split into sums over different flavours Then, we need to map the emitter-spectator pairs of the virtual sector into the dipoles of the real sector. Each pair (i, j) can be mapped in more than a single dipole (a, c, b), where parton ac is the emitter, c is the radiation and b is the spectator. Therefore, to distinguish the case in which the pair (i, j) is linked to a specific dipole (a, c, b), we use the following notation for the momentum mappings of Section 3.1 The application of Eq.(A.4) in the corresponding dual counterterm leads to We now want to use Eq.(A.5) and count the number of counterterms through the number of dipoles in the real sector, by transforming each sum on the right-hand side of Eq.(A.3) in the following way where #(f ) m denotes the number of partons with flavour f in the m-parton configuration and #{a, c} m+1 the number of pairs {a, c} in the (m + 1)-parton configuration such that the emitter parton ac has flavour f . Given an m-parton configuration, all the possible (m + 1)-parton processes can be obtained either by increasing by one the number m g of gluons, or by decreasing by one m g and adding a quark-antiquark pair. Let us focus on the first case, which includes all the dual counterterms beside the quark contributions to the gluon wave-function renormalization. The latters are required in the second class of (m + 1)-parton processes, treated at the end of this Appendix. Therefore, for the moment we assume that the counterterms C g,j only include gluon and ghost contribution to the gluon renormalization. We have m q and mq being the number of quark and antiquarks of a given flavour, respectively. The counting in Eq.(A.6) is possible only if the integral of C i,j over the m-particle phase space does not change when we exchange parton i with a parton of the same flavour. This condition is certainly fulfilled if we use the q 1 -cut method and so always set on-shell the loop momentum propagator connecting the external partons (q 1 ). However, particular attention must be payed to this point if we want to use the q 1 -q 2 -method. In fact, if a and a had the same flavour and we used σ ac,b to describe the radiation collinear to a and σ D (2) bc,a for the emission collinear to a , then the dual subtractions would differ among each other.
As a solution to this problem we can perform, for each dual subtraction, the following symmetrization σ D (1) which ensures that the dual subtractions depend only on the flavours of the partons in the associated dipoles. Eq.(A.8) may be derived from a different prescription for the application of LTD, which can be summarized as follows: when we consider the generic contribution A (be it a triangle C 0 , a bubble B 0 or a wave-function renormalization counterterm ∆Z) from a given emitter-spectator pair (i, j) of the virtual sector, we first use the identity A = A/2 + A/2; then, for the first term, we apply LTD as usual while, for the second term, we denote the internal momenta e we apply LTD as if we had considered the pair (j, i) instead of (i, j). From a diagrammatic point of view, this means that we associate the triangle in Fig.(3a) to the first A term and the triangle in The presence of a factor 2 in front of the sum over gluon-gluon pairs on the right-hand side of Eq.(A.10) has a clear interpretation: given a certain pair of gluons in the real sector, their roles can be exchanged and, therefore, we need to take twice the dual counterterm. In this way, indeed, in one counterterm we can apply the mapping where one of the two gluons is the radiation and, in the other counterterm, we can use the mapping where the role of radiation is played by the other gluon. Furthermore, using the relation The sums in Eq.(A.12) present one dual counterterm for each dipole of the real phase space. Therefore, Eq.(A.12) tells us how to associate every single counterterm of the dual subtraction scheme with the proper singular configuration of the real phase space. Let us now move to processes with a quark-antiquark pair replacing a gluon. These involve the counterterms C g,j , that we assume to include only the quark contribution to the gluon renormalization, since gluon and ghost contributions are used in processes with one additional gluon in the final state. Before to apply Eq.(A.6), we first multiply and divide the counterterms by N f and then use the N f in the numerator to introduce a sum over the different quark flavours i=g j =i

B Dual Counterterms for the massless case
In this section we collect all the formulas to be used for a numerical implementation of the dual subtractions. To regulate the divergences of the real matrix elements, we have cut the relevant loop integrals and collected their divergences. Then, in principle, we could organize the integrated counterpart in the form of finite remnants of the loop integrals. This could be useful for implementations of our formalism in general purposed programs that perform the full reduction of the virtual amplitude. On the other hand, to make the implementation easier, in this appendix we will continue to treat our cut loop integrals as normal countertems and proceed to present the list of their integrated version. The formulation of our method is given as follows: in Table 1 we list the equation number of the dual subtraction along with the name of the corresponding jacobian factor that converts the radiation variables to that of the loop momentum as explained in Section 3. We label with p a and p b the partons in the initial state and p i , p j , p k and so on for the ones in the final state. Moreover, we remember here that the formulae for the dual dipoles are written following the convention of in-going and out-going momenta, p a + p b = p i + p j + p k + ... . Furthermore, the names of the counterterms are written using the same conventions as in the presentation of the Altarelli-Parisi splitting functions in Section 3.2. The prefixes specify the emitter and spectator respectively. The suffixes have to be understood as (ac)c for a final state emitter (ac) radiating the parton c and it is a(ac) for an initial state emitter (ac) (the hard parton) that is the result of the splitting of the incoming parent parton a. So for example: if _qq stands for the counterterm for an initial state quark (i in the first part of the prefix and q in the first part of the suffix) emitting a gluon and becoming the hard quark of the underlying event (the second part of the suffix q), and this gluon is connected to a final state spectator (as indicated by the f in the second part of the prefix); f f _qg stands instead for a virtual quark (q) splitting into a qg pair with the gluon (g) being the radiated parton connected to another final state parton. The dual dipoles formulae do not depend from the location of the spectator momentum. The formulae for the dipole subtractions do not depend from the specific momentum mappings nevertheless for their numerical implementation one choice has to be made. Here we adopt the momentum mappings reported in ref. [5] that from now on we dub CS and will refer to the equation x.y in that paper with CS(x.y) for easy of reference 4 .

B.1 final-final
For the case of final state emitter and spectator we will report in detail all the formulas that we take from CS. For the other cases we will just point to the other formulas in the same paper. The momentum mapping to transform the three momentum {i, j, k} into {ĩj,k} is defined in the Equation The factorization of the phase space in Eq. CS (5.17-20) is then dφ(p i , p j , p k ; Q) = dφ(p ij ,p k ; Q)[dp j (p ij ,p k )] (B .3) where [dp j (p ij ,p k )] = d d p j (2 π) d−1 δ + (p [dp j (p ij ,p k )] = (2p ijpk ) 1−ε 16 π 2 dΩ d−3 (2 π) 1−2ε dz i dyΘ(z i (1 − z i ))Θ(y(1 − y)) × (z i (1 − z i )) −ε (1 − y) 1−2 ε y −ε . (B.6) As explained in Section 3, the local cancellation of the singularities among virtual and real contribution is obtained only if one takes into account also the jacobian factor associated to the mapping. In particular we have to divide the formulae for the counterterms in Tab.1 by the factor I(p j ;p ij ,p k ) that for this case amount to multiply them by Furthermore, instead of the whole region 0 < z i , y < 1 implied by the Θ functions in Eq.(B.6) the restriction of the phase space induced by the R 1 function corresponds to limit the integration region to the configurations with p i p j < p j p k or in terms of the radiation variables to the restriction We introduce the normalization factor N ε = (4 π) ε Γ(1 − ε) (B.9) and the short hand notation l 2 ≡ log(2). We integrate the counterterms built by the product of the dual integrands multiplied times the corresponding jacobian factor (jac(ff) in this case), using the measure in Eq.(B.6) and restricting the integration to the region in Eq.(B.8), obtaining Note that we adopt the convention to have in the prefactor the kinematic invariant built with the momenta of the emitter and spectator partons. The same convention is adopted also for the other cases [72]. Following the approach of combining the integrated countertems with the virtual matrix element it is more convenient to do so because this is the invariant which is kept fixed during the x integration. In the above equations the case for ff_gq_int is included in ff_gq_int as a factor of 2 compensated with the T R factor. Furthermore, the two parts corresponding to the simmetrization in Eq.(3.39) gives identical integrated contribution. We point out that the momentum mapping in Eqs.(B.1,B.2) is symmetric for the exchange of p i and p j and so there is just one counter event although its weight is determined by the phase space restrictions that event by event are different for the two contributions. The same considerations also apply to the next case that is the final-initial.

B.2 final-initial
The momentum mapping to transform the three momentum {i, j, a} into {ĩj,ã} is defined in the Equation CS(5.37)p in terms of the variables x = p i p a + p j p a − p j p i p i p a + p j p a (B.24) u j = p j p a p j p a + p i p a = z j .
This momentum mapping coincides with the previous one. This is particularly useful for our approach because such a map connects a set of radiation variables and a Born configuration with an initial(final) state emitter and a final(initial) state spectator to a unique real configuration and vivecersa. The factorization of the three parton phase space in Eq. CS(5.70-73) is the same as in the previous case and the jacobian factor needed to transform the radiated momentum into the loop momentum is again given by jac(if) = 1 − u j = z i = jac(fi) (B.25) The integration region corresponding to p a p j < p i p j in terms of the radiation variables is now given by We integrate the counterterms built by the product of the dual integrands multiplied times the corresponding jacobian factor in Eq.(B.25), using the measure in Eq.(B.17) and restrict-