Four-dimensional unsubtraction with massive particles

We extend the four-dimensional unsubtraction method, which is based on the loop-tree duality (LTD), to deal with processes involving heavy particles. The method allows to perform the summation over degenerate IR configurations directly at integrand level in such a way that NLO corrections can be implemented directly in four space-time dimensions. We define a general momentum mapping between the real and virtual kinematics that accounts properly for the quasi-collinear configurations, and leads to an smooth massless limit. We illustrate the method first with an scalar toy example, and then analyse the case of the decay of a scalar or vector boson into a pair of massive quarks. The results presented in this paper are suitable for the application of the method to any multipartonic process.


Introduction
The development of new computational techniques to obtain more accurate theoretical predictions at colliders has been strongly pushed forward by the high precision experimental data obtained from the LHC. In the framework of perturbative quantum field theory (and perturbative QCD in particular), the presence of infrared (IR) and ultraviolet (UV) divergences constitutes the main difficulty that must be overcome to obtain physical results. renormalisation has been proven to deal successfully with the UV structure of these theories, and the cancellation of IR singularities is guaranteed by general theorems [2,3] for physical observables that sum over all the possible degenerate IR configurations. This requires taking into account both loop scattering amplitudes and real processes with the emission of additional external particles; the sum of all contributions leads to IR-safe observables. On the other hand, UV divergences are removed by suitable counter-terms, whose divergent structure depends only on the specific theory under consideration.
In order to make these divergences manifest explicitly, the standard approach relies in the introduction of a convenient regularisation method. Nowadays, the standard choice in gauge theories is dimensional-regularisation (DREG) [4,5,6,7] because it preserves gauge invariance. Within DREG, the space-time is analytically continued from d = 4 to d = 4 − 2 dimensions; both UV and IR divergences appear as -poles. Using as a regulator, it is possible to perform the loop and the phase-space integrals for the virtual and real-radiation amplitudes, respectively. Thus, the poles of the virtual corrections are cancelled with those present in the real-radiation contributions (due to soft/collinear configurations) and those included in the UV counter-terms. The usual procedure in this framework is to regularise each contribution separately, integrate the expressions and, finally, cancel the poles and take the limit → 0.
Besides the renormalisation and the regularisation method, the fact that real and virtual contributions have the same IR-divergent behaviour is the underlying basis of the subtraction methods [8,9,10,11]. The general idea of subtraction is the introduction of counter-terms which mimic the local IR behaviour of the real components and that can easily be integrated analytically. In this way, the integrated form is combined with the virtual component whilst the unintegrated counter-term cancels the IR poles originated from the phase-space integration of the real-radiation part. The subtraction paradigm has evolved to different versions from the FKS-subtraction [8,9], and dipole-subtraction [10,11], to antenna-subtraction [12], q T -subtraction [14,15] and other recent variations [16,17,18,19,20,21,22]. The treatment of massive particles has also been considered specifically [23,24,25,26,27]. However, all these approaches treat separately real and virtual corrections, since the final-state phase-space of the different contributions involves different numbers of particles. The construction of IR counter-terms is inherent to the subtraction approach, and it constitutes the main restriction for an efficient application to multi-leg multi-loop processes.
With the purpose of obtaining a fully-local cancellation of IR singularities without building IRcounter-terms, we apply the loop-tree duality (LTD) theorem [28,29,30,31] to manage the virtual corrections. This theorem asserts that loop integrals are expressible as the sum of dual integrals, which are built starting from tree-level like objects and replacing the loop measure with an extended phasespace measure. The main advantage of LTD is that by introducing additional physical on-shell particles, dual integrals and real-radiation contributions exhibit a similar structure. Moreover, the loop threshold and IR singularities are always restricted to a compact region of the loop three-momentum space [32,33]. These two properties of LTD allow a natural integrand-level combination of virtual and real corrections. The method has been recently developed for processes involving massless particles in Refs. [36, 37, 38, are labeled as p i with i ∈ {1, 2, . . . N }, then the internal virtual momenta are given by with the loop internal momentum and k N = 0 due to momentum conservation. The corresponding expression for the scalar integral is where is the scalar Feynman propagator associated to a virtual particle with mass m i and four-momentum q i,µ = (q i,0 , q i ) (q i,0 is the energy and q i are the spatial components). We recall that the +ı0 prescription is introduced to separate, in the imaginary axis, the solutions arising from the on-shell condition, i.e. G F (q i ) −1 = 0. In particular, this translates into two solutions with positive and negative imaginary components, respectively. On the other hand, denotes the standard one-loop integration measure in d-dimensions. According to the LTD theorem, any loop contribution to scattering amplitudes in any relativistic, local and unitary quantum field theory can be computed through dual integrals, which are built from single cuts of the virtual diagrams at one-loop [28]. In other words, there exists a formal connection among loop and phase-space integrals. The cut condition is implemented by restricting the integration measure through the introduction ofδ which forces to integrate in the positive energy mode (q i,0 > 0) of the corresponding on-shell hyperboloid (i.e. q 2 i = m 2 i ). It is worth appreciating that, in the massless limit, all the hyperboloids associated with the on-shell condition G F (q i ) −1 = 0 degenerate into light-cones. Considering the one-loop scalar integral, the LTD theorem establishes that its dual representation is given by i.e. the sum of N dual integrals, each one associated with one of the possible single-cuts. In Eq. (6), the dual propagators are with i, j ∈ {1, 2, . . . N }, k ji = q j − q i and η an arbitrary future-like or light-like vector, η 2 ≥ 0, with positive definite energy η 0 > 0 * . Since η is arbitrary, we can chose η µ = (1, 0) to simplify the computations.
Assuming that there are only single powers of the Feynman propagators, the dual representation in Eq. (6) is straightforwardly valid for loop scattering amplitudes. The single-cuts do not affect numerators, therefore, the dual representation of scattering amplitudes is obtained by simply adding all possible dual single-cuts of the original loop diagram, and replacing the uncut Feynman propagators by dual propagators. If there are higher-powers of the propagators, however, we should apply the extended version of the LTD theorem [31] by using the Cauchy's residue theorem with the well-known formula for poles of order n, i.e.

Res(A, q
where q (+) i,0 = q 2 i + m 2 i − ı0 is the positive energy solution of the corresponding on-shell dispersion relation. In that case, the explicit form of the scattering amplitude is relevant because the numerator is affected by the derivative.

Massive scalar three-point function within LTD
We present in this Section the first analytical application of LTD with massive particles. In particular, we consider the scalar three-point function with one massless internal state and the remaining internal and outgoing particles with mass equal to M . The final-state on-shell momenta are labeled as p 1 and p 2 , with p 2 1 = M 2 = p 2 2 , and the incoming one is p 3 = p 1 + p 2 ≡ p 12 , by momentum conservation, with virtuality p 2 3 = s 12 . We consider s 12 > 0, i.e. we work in the time-like (TL) kinematical region. The internal momenta are q 1 = + p 1 , q 2 = + p 12 and q 3 = , where is the loop momentum (see Fig. 3). When the internal lines are set on-shell, the momentum q 1 is massless, whilst the other two are massive and fulfill q 2 2 = M 2 = q 2 3 . This is the master scalar configuration for the calculation of the QCD corrections to the physical case A * → qq(g) with massive quarks that will be studied in Sec. 9. We define as the normalised mass and velocity, respectively. The well-know result of this massive scalar three-point function is given by [52,53] L (1) with and c Γ the one-loop volume factor Applying LTD, the dual representation of the scalar integral in Eq. (10) consists of three contributions with The corresponding on-shell hyperboloids are shown in Fig. 2 (left). Due to the rotational symmetry, it is enough to show the ( 0 , z ) plane. As discussed in Ref. [32], the intersection of on-shell hyperboloids is associated with multiple internal states becoming simultaneously on-shell. In this case, the forward on-shell hyperboloids of G F (q 1 ) and G F (q 2 ), and the backward one of G F (q 3 ) intersect in a single point: this leads to a soft singularity. The other intersection takes place between the forward on-shell hyperboloid of G F (q 2 ) with the backward of G F (q 3 ), which manifests as a threshold singularity in I 2 .
Notice that there are not collinear singularities, because the mass prevent that the on-shell hyperboloids degenerate into light-cones. In that situation, there would be extended forward-backward intersections in the ( 0 , z ) plane leading to collinear poles, as described in the massless case studied in Refs. [36,39]. Now, we shall explicitly compute these dual integrals. We choose first a proper reference frame to simplify the analytic expressions. Hence, we work in the centre-of-mass frame of p 1 and p 2 , and parametrise the involved momenta as Figure 2: On-shell hyperboloids of the massive three-point function in the loop coordinates µ = √ s 12 /2 (ξ 0 , ξ x , ξ y , ξ z ) in two dimensions (left plot); forward and backward on-shell hyperboloids are represented by solid and dashed lines, respectively. The intersection of on-shell hyperboloids leads to soft and threshold singularities in the loop three-momentum space (right plot), collinear singularities are regulated by the mass.
In order to describe the internal on-shell momenta, we must take into account that q 1 corresponds to a massless state (q 2 1 = 0), whilst q 2 and q 3 are associated with massive virtual particles (q 2 2 = M 2 = q 2 3 ). For this reason, we write where ξ 1,0 , ξ 2 , ξ 3 ∈ [0, ∞) and v i ∈ [0, 1] are the integration variables describing the modulus of the three-momentum and polar angle of the loop momenta, respectively. With these variables, the LTD transforms the loop integration measure into for the massive case (i = 2, 3), whilst it reduces to for a massless state (i = 1). Notice that ξ 1,0 = ξ 1 since q 1 is massless. The integration of the loop momentum in the transverse plane, which is described by the unit vectors e i,⊥ is trivial. The scalar products of internal with external momenta are given by which reduce to 2q 1 · p 1 /s 12 = ξ 1,0 v 1 and 2q 1 · p 2 /s 12 = ξ 1,0 (1 − v 1 ) for a massless on-shell state. Using these variables, the dual integrals in Eq. (14) are rewritten as , .
Notice that I 2 contains a threshold singularity at m 2 + ξ 2 2 = 1, i.e. ξ 2 = β ≤ 1. These integrals can be calculated analytically to all orders in in the massless limit [36]. In this limit, they read with c Γ = c Γ / cos(π ) the phase-space volume factor. As expected, the sum of the three dual integrals in Eq. (22) agrees with the well-known massless scalar three-point function.
The massive case is a bit more cumbersome because of the dependence on m. Again, I 1 vanishes because the energy integral factorises and it lacks of any characteristic scale. Actually, I 1 is singular both in the IR and UV. However, the sum of the three dual integrals and the equivalent original Feynman integral contain only soft divergences. The other two dual integrals can be integrated in the angular variable analytically, thus keeping the exact -dependence. This leads to However, an expansion in is necessary to integrate in the modulus of the loop three-momentum, which leads to a final result that includes corrections up to O( 0 ). Besides that, the two dual integrals Eq. (23) are singular only in the UV. Therefore, we introduce the following expansion with g UV = lim ξ i →∞ g i (ξ i ). The first term in the r.h.s. of Eq. (24) gives the same result for both dual integrals, i.e.
For the second term, which is regular in the UV, we perform the change of variable with z ∈ [1, ∞). The total result for the real part of the two dual integrals, up to O( ), reads The dual integral I 3 is purely real, while the dual integral I 2 generates an imaginary component due to the intersection of the forward on-shell hyperboloid of G F (q 2 ) with the backward on-shell hyperboloid of G F (q 3 ). Its imaginary part can be calculated to all orders in from which is the expected result obtained through the application of the Cutkosky's rule. The sum of these contributions in Eq. (27) and Eq. (28) agrees with the original Feynman integral in Eq. (10).

Massive scalar decay rate in DREG
In order to establish a physical parallelism and understand the subtraction of IR singularities, we work in a simplified toy scalar model with a massive scalar particle φ which couples to a massless one, ψ. In concrete, we consider the decay process φ(p 3 ) → φ(p 1 ) + φ(p 2 ), with p 2 1 = M 2 = p 2 2 (on-shell massive external particles) and p 2 3 = s 12 (off-shell incoming particle). The Born-level decay rate is given by where g is the coupling and s 12 > 4M 2 to guarantee the physical feasibility of the process. To compute the corresponding NLO correction, we need to add virtual (i.e. one-loop) and real (i.e. extra-radiation) contributions. We will assume the presence of only one massless particle inside the loop, as well as the emission of a massless real particle in the extra-radiation contribution. The corresponding NLO diagrams are exhibited in Fig. 3 † .
Let's start with the virtual part, which we assume proportional to the scalar three-point function, i.e. Figure 3: Kinematic configuration of the NLO corrections to the decay process φ → φ+φ. The one-loop contribution is proportional to the scalar three-point function, with a virtual massless ψ inside the loop (left). The real contribution is due to the interference terms originated by the emission of an on-shell massless particle ψ (right). In this case, the momentum configuration is given by Since s 12 > 0, the virtual decay rate is given from Eq.
In order to calculate the totally inclusive decay rate, we have to consider the real emission process. In this particular toy-example the cancellation of IR-singularities is achieved by including only the interference terms originated in the process φ(p 3 ) → φ(p 1 ) + φ(p 2 ) + ψ(p r ), as depicted schematically in Fig. 3. Explicitly, where we used the definition of the massive three-body phase space in Eq. (100) and y ir = 2 p i · p r /s 12 .
To compute this integral, we apply the change of variables suggested in Appendix A, which allows to factorise the energy and the angular dependence of the integrand. By using Eq. (102), we obtain The integration in w can be trivially performed, and it leads to the appearance of an -pole. The integral in z is finite if x S > 0 (i.e. in the massive case), thus we expand the integrand in before the integration. The resulting expression is Putting together the virtual and real contributions from Eq. (31) and Eq. (34), we get with a = g 2 /(4π) 2 . The purpose of the following discussion will be the derivation of a purely fourdimensional representation of this result, through the local cancellation of all the IR divergences present in both real and virtual contributions. It is crucial to emphasise that this cancellation of IR singularities at integrand level is achieved by a suitable mapping of momenta.
5 Phase-space partition and real-virtual mapping with massive particles The first ingredient that we need to introduce is a complete partition of the real phase-space [36,39], in such a way that each individual region of that partition contains a single soft, collinear or quasicollinear configuration. The quasi-collinear configurations are those in which a massless particle becomes collinear with a massive one [23]. In that case the mass acts as an IR-regulator and the collinear -poles that appear in DREG in the massless case are transformed into finite logarithmic terms in the mass. These logarithmic contributions are cancelled in the total cross-section but the massless and the → 0 limits do not commute for the virtual and real corrections separately. Thus, we split the real phase-space by defining the domains For instance, R i selects configurations with p i p r or close to collinear, excluding the remaining ones.
In particular, this partition reduces to for the simplest 1 → 3 scenario. The definition of the phase-space partition in Eq. (36) is the same that we would use in the massless case [39]; the only difference is that now some of the external momenta are massive. That partition, together with the well-motivated mapping of momenta that will be presented in the following, ensures a smooth massless limit, therefore an integrand level cancellation of the logarithmic dependences in the mass arising from the quasi-collinear configurations of the virtual and real corrections, and thus a more stable numerical implementation of the method.
Then, we shall define a proper momentum mapping in each region to match the singular behaviour of the real and the dual integrands. In order to properly combine real and virtual contributions at integrand level, we need to generate the N + 1 on-shell kinematics by making use of the N -parton Born-level process and the on-shell loop momenta. One of the main difficulties in constructing a momentum mapping with massive particles is that the on-shell conditions lead to quadratic equations in the mapping parameters. However, it is very well-known that massive vectors can be expressed in terms of two massless momenta. Then, we can exploit this property to simplify the mapping equations. For the case of a pair of particles of the same mass, the corresponding massive momenta can be written as withp 2 1 =p 2 2 = 0 and β ± = (1 ± β)/2. Moreover, the massless momenta fulfil the following useful identities In their centre-of-mass frame, these massless momenta are simply given bŷ Going back to the toy example in Sec. 4, let's start with the first region where R 1 = 1 (i.e. y 1r < y 2r ). Motivated by the factorisation properties of QCD in the quasi-collinear limit and the momentum decomposition in Eq. (38), we propose the following mapping with q 2 which fulfils the momentum conservation constraint by construction. As in the massless case, the momentum p 2 acts as the spectator of the splitting process and is used to balance momentum conservation. The emitters have momenta p 1 and p 1 , and have the same mass. Although restricted to three final-state particles, the momentum mapping in Eq. (41) can easily be generalised to the multipartonic case, with p k = p k for k = 1, 2, r. The parameters α 1 and γ 1 are determined from the two on-shell conditions whose explicit solutions are whilst (p r ) 2 = 0 by construction, since q 2 1 = 0. Due to the fact that we are dealing with quadratic equations, there are two sets of solutions. The solution in Eq. (43) is compatible with the soft limit; it recovers the Born-level kinematics when ξ 1,0 → 0. In that limit, (α 1 , γ 1 ) → (β − , β + ) and therefore (p 1 , p 2 ) → (p 1 , p 2 ). Also, it properly reduces to the massless parametrisation defined in Refs. [36,39], i.e. if m → 0, we have Using these definitions, the kinematical invariants y ij become Figure 4: The dual integration regions in the loop three-momentum space, with ξ ⊥ = ξ 2 x + ξ 2 y .
which fullfil y 12 + y 1r + y 2r = 1 − m 2 /2. Again, we recover easily the massless expressions [36,39] with α 1 = 0. In order to improve the presentation of the results, it is also convenient to express the mass in terms of α 1 , Then, we compute the associated Jacobian in the physically allowed region (i.e. those points belonging to the domain R 1 ), which is given by with dy 1r dy 2r = J 1 (ξ 1,0 , v 1 ) dξ 1,0 dv 1 . Notice that this expression is apparently free of square roots, since the mass dependence was rewritten in terms of α 1 , as suggested in Eq. (46).
On the other hand, we have to express R 1 in terms of the dual variables. If we use the mapping given in Eq. (41), we obtain which is the characteristic function associated to the domain R 1 . In Fig. 4, we show this domain in the loop three-momentum space. The massless limit agrees with the expected result. Moreover, the three-body phase-space limits defined by the condition h p = 0 are simply determined by v 1 = 0.
In the complementary region, where R 2 = 1 (i.e. y 2r < y 1r ) with q 2 2 = M 2 , the mapping is defined by ‡ where are the associated on-shell conditions. In this case, the condition (p 2 ) 2 = M 2 is fulfilled by construction.
Solving the system and selecting the physical solution, we get In order to check the consistency of this solution, we consider the massless limit, obtaining which implies that the parametrisation reduces to the expected one. On the other hand, the two-body invariants are given by and, since this mapping will be used in the region of the real phase-space defined by R 2 , we rewrite the associated characteristic function as The corresponding domain in the loop three-momentum space is also shown in Fig. 4. The associated Jacobian is given by with dy 1r dy 2r = J 2 (ξ 2 , v 2 ) dξ 2 dv 2 , and we made use of the identity to simplify the expressions.

General momentum mapping
The momentum mappings previously presented can easily be extended to the most general multipartonic case in which the emitter and the spectator have different masses: p 2 i = m 2 i and p 2 j = m 2 j , respectively. The decomposition of their momenta in terms of two massless momenta (p 2 i =p 2 j = 0) is given by with where is the usual Kallén function. The massless momenta fulfil the useful conditionp i +p j = p i + p j . The mapping with the momenta of the real process is formally equal to the mapping already considered in Eq. (41), i.e.
It leads to the on-shell conditions In Eq. (60), we have imposed that the spectator and the radiated particle have the same flavour (and, thus, the same mass) in the virtual and real processes; p 2 j = (p j ) 2 = m 2 j and q 2 i = (p r ) 2 = m 2 r , respectively. The emitter, however, might change flavour, (p i ) 2 = (m i ) 2 = m 2 i . This situation occurs, for instance, when a gluon splits into a massive quark-antiquark pair. The solution to Eq. (60) for the parameters of the mapping reads with The momentum mapping in Eq. (59) has an smooth limit whenever any of the involved particles becomes massless; in particular, if the spectator is a massless particle, then α i = 0.

Massive scalar decay rate from four-dimensional unsubstraction
In this section we illustrate the method of LTD four-dimensional unsubstraction [36,39] with the massive toy example presented in Sec. 4. All the necessary ingredients have been presented in the previous sections. We combine at integrand level the dual loop contributions (Sec. 3) with the real-radiation terms (Sec. 4) with the help of the momentum mappings defined in Sec. 5. Since the sum of all the contributions is UV and IR finite, the final result is free of -poles. We would like to emphasise that, in a generic situation, this assertion is not enough to guarantee the integrability of the expressions in fourdimensions. However, by virtue of the momentum mappings and the unification of the dual coordinates, LTD leads naturally to a local cancellation of divergences and the limit → 0 can be considered at integrand level.
The LTD representation of the virtual decay rate in this toy example is given by with where the dual integrals I i are defined in Eq. (21), and we take the integration measure exactly with = 0. In order to ensure the cross-cancellation of spurious singularities and get a direct = 0 limit we must rewrite all the on-shell momenta in terms of the same coordinate system. This change of variables is explained in Appendix B. With that change of variables the virtual decay rate in Eq. (63) becomes a single unconstrained integral in the loop three-momentum.
We also consider the real contribution given by Eq. (32). First, we split the real three-body phasespace according to Eq. (37), and define that obviously fulfill Γ Second, in each region of the real phase-space we apply one of the momentum mappings defined in Eq. (41) and Eq. (49), respectively. The main advantage of these mappings is that they are optimised to deal smoothly with the massless limit in each of the two regions. Thus, we rewrite the real contributions in terms of the loop variables and we obtain where the Jacobians of the respective transformations are given by Eq. (47) and Eq. (55), and the integration domains, which are restricted by the functions R i from Eq. (48) and Eq. (54), are shown in Fig. 4. Again, after applying the change of variables defined in Appendix B to bring the two real contributions to a common coordinate system, we can consistently take the limit = 0, directly at integrand level.
The sum of all the virtual and real contributions, Eqs. (63), (67) and (68), is a finite function in the = 0 limit because all the IR singularities are cancelled locally in the loop three-momentum space at integrand level. The virtual contribution, however, contains a threshold singularity at ξ 2 = β. This singularity is integrable and can be treated numerically by contour deformation [33,34,35]. For the toy scalar model and the physical examples that we are considering in this article there is a simplest solution: we can compactify the high-energy region with ξ > β into the unit sphere by using a change of variables. Explicitly, since the integrand is a function of the modulus of the three-momentum and the polar angle, then where the threshold singularity has been mapped into to the upper end-point, namely x = 1. This approach is very efficient for the numerical implementation.
Finally, we numerically integrate simultaneously the virtual and real corrections from Eqs. (63), (67) and (68) with the help of Eq. (69) to obtain the total decay rate at NLO, Γ (1) , as a function of the dimensionless mass parameter m. The result is shown in Fig. 5, and it is compared with the DREG analytic expression given by Eq. (35). The agreement is excellent and quite stable numerically. Moreover, the massless transition is very smooth because the momentum mappings are optimised to deal with the quasi-collinear configurations. In other words, the massless limit can directly be taken at the integrand level. This is another interesting advantage of the LTD approach.

Unintregrated wave function and mass renormalisation for heavy quarks
In order to consider physical processes with heavy quarks we should also take into account self-energy corrections. The well-known expressions of the wave function and mass renormalisation constants, in the Feynman gauge with on-shell renormalisation conditions, are not suitable, in particular, for the implementation of a local subtraction of the IR singularities. We shall provide unintegrated expressions. The case of massless quarks has been studied in detail in Ref. [39]. In Eq. (70), we explicitly identify the origin of the -poles, and the IR singularities of the wave function should cancel the IR singularities arising from the squared amplitudes of the real processes with radiated gluons off quarks.
We consider the process in which there are two on-shell massive fermions with momenta p 1 (quark) and p 2 (antiquark). The explicit one-loop self-energies are given by In these expressions, we keep the same internal momenta q i that were used to define the vertex corrections in Sec. 3 (See also Fig. 8 in Appendix C). This will allow us to reuse the same momentum mappings already defined to treat the vertex corrections. Because of the symmetry p / 1 ↔ −p / 2 , the unintegrated expression for the antiquark self-energy corrections can be deduced from those of the quark. Thus, we consider in the next only the quark self-energy. In the Feynman gauge Working in the on-shell renormalisation scheme (OS), the renormalised self-energy fulfills from where the wave function and mass renormalisation corrections are given by Explicitly, from Eq. (73) we obtain It is worth to stress that the expression of the wave function renormalisation constant in Eq. (76) tends smoothly in the massless limit to the corresponding expression given in Ref. [39]. It is also relevant to notice that the term proportional to M 2 (G F (q 3 )) 2 leads to soft divergences when q 1 gets on-shell. They are expected to cancel the soft divergences of the squared amplitudes of the real corrections. The dual representation of the mass renormalisation constant is straightforward from the LTD theorem. The term M 2 (G F (q 3 )) 2 in Eq. (76), however, introduces double poles that need to be treated specifically [31]. The final dual representations for both renormalisation factors are or in term of dual variables defined in Sec. 3

UV renormalisation
We shall now remove the UV divergences of the renormalisation constants by defining suitable integrand level UV counter-terms. These UV counter-terms are obtained by expanding Eqs. (76) and (77) around the UV propagator G F (q UV ) = 1/(q 2 UV −µ 2 UV +i0), where q UV = +k UV with k UV arbitrary [45,36,39].
The simplest choice is k UV = 0. We obtain whose integrated form is The sub-leading terms in Eq. (80), which are proportional to µ 2 UV , have been adjusted in such a way that only the UV poles in Eq. (70) are subtracted at O( 0 ). Therefore only contain IR singularities, including the finite terms which are scheme dependent.
The dual representation of Eq. (80) requires to evaluate the residue of poles of second and third order [31,36,39] located at q (+) UV,0 = q 2 UV + µ 2 UV − ı0. We obtain Then, we use the parametrisation with m UV = 2µ UV / √ s 12 , and the UV counter-term get the form , .
Similarly, we should subtract the UV singularities of the qqA interaction vertex, with A = {φ, γ, Z} for the explicit examples that we will consider later. As for the self-energy contributions, the UV counterterm is obtained by expanding the vertex corrections around the UV propagator G F (q UV ) = 1/(q 2 UV − µ 2 UV + i0). In the Feynman gauge, the generic expression of the vertex UV counter-term reads where the the tree-level vertices Γ A are given in Eq. (109) of Appendix C. The µ 2 UV term is sub-leading and the coefficient d UV is adjusted to subtract only the UV pole. Performing the explicit calculation, we find that in the MS scheme these coefficients are Notice that this choice of the sub-leading contributions differs from d γ,UV = 4 proposed in Ref. [45]. The difference is, however, of O( ). The integration of the vertex UV counter-term leads to the result with that translates into The dual representation of Eq. (86) is (see Ref. [39]) After an explicit calculation, we obtain for the vertex UV counter-terms where the function integrates to zero and does not contribute to the renormalisation of the vertex. However, this additional term is necessary to achieve a local cancellation of the UV behaviour.
The UV divergences of the wave function cancel exactly the UV divergences of the vertex corrections for photons and Z bosons, because conserved currents or partially conserved currents, as the vector and axial ones, do not get renormalised. The corresponding dual representations, however, do not cancel each other at integrand level. In particular, the wave function renormalisation contains linear UV singularities that cancel upon integration. Also, the vertex UV counter-term contains terms that are proportional to the mass and cancel upon integration. The contribution of all this spurious terms is, however, crucial to cancel locally all the UV singularities. The coupling to scalar particles, on the contrary, needs to be renormalised

LTD four-dimensional unsubtraction for physical processes
We have already defined all the necessary ingredients to test the four dimensional implementation of NLO corrections to physical processes in the LTD framework. In particular, we will compute the NLO QCD corrections to the decay rate A * → qq(g), with A = φ, γ, Z. The actual implementation is indeed independent of the decaying particle. The renormalised one-loop amplitude is given by where |M is the unintegrated UV counter-term of the one-loop vertex correction, |M A , and Z IR 2 (p i ) are the IR components of the quark and antiquark self-energy corrections. From the renormalised one-loop amplitude |M , which contains only IR singularities, we construct the LTD representation of the renormalised virtual decay rate The corresponding dual amplitudes for the vertex corrections are given explicitly in Eqs. (115), (116) and (117). As for the toy scalar example presented in Sec. 6, the real contributions are implemented by splitting the real phase-space in two domains with Γ (1) R,A,2 the real total decay rate. The real emission squared amplitudes are given in Eq. (113). In each of the real-phase space domain we introduce one of the momentum mappings defined in Sec. 5. The sum of the virtual and real corrections in Eq. (96) and Eq. (97) is a single integral in the loop three-momentum. It is UV and IR finite locally and thus can be calculated numerically with = 0. We follow the same numerical implementation as for the scalar example presented in Sec. 6, and compare the numerical output with the analytical total decay rate, which has the form Our results, normalised to the LO decay rate Γ A , are presented in Fig. 6 for a scalar and pseudoscalar, and in Fig. 7 for vector bosons. The agreement with the analytic prediction is excellent in all the cases. Moreover, the massless limit, i.e. x S → 0, is also well-defined and we recover the known results. This is a very subtle point, because individual contributions in DREG are not smoothly well-defined in that limit.

Conclusions and outlook
In this article, we have generalised the four-dimensional unsubtraction method [36,37,38,39,40] to deal with massive particles. Based in the LTD theorem, it exploits the possibility of expressing virtual amplitudes as phase-space integrals, and the fact that threshold and IR singularities are always restricted to a compact region of the loop three-momentum integration domain. This is a crucial point, because it allows to establish a momentum mapping to generate the real-emission on-shell kinematics starting from the Born level momenta and the loop three-momentum, in such a way that a local cancellation of both IR and UV singularities can be achieved without introducing IR subtractions. In particular for the massive case, we have defined a general momentum mapping that accounts properly for the quasicollinear configurations.
First, we started by inquiring in the computation of the scalar three-point function with massive particles within the LTD approach. Besides recovering the previously known results, the analysis of the integration domains of the dual contributions allowed us to understand the origin of its singular structure. Then, we illustrated the local cancellation of IR and quasi-collinear configurations with a toy scalar example. /Γ (0) Figure 6: Total decay rate at NLO for scalar and pseudoscalar particles into a pair of heavy quarks as a function of the mass, normalised to the LO. In the left panel, we consider a standard Higgs boson, whilst in the right panel we plot the decay rate for a pseudoscalar particle (c q = 1/2). The solid blue lines correspond to the usual DREG analytic result, while the red dots were computed numerically within the LTD unsubtracted method. We also consider a renormalisation scale variation, in the range The full cancellation of IR singularities requires the contributions of the self-energy corrections. Thus, we defined unintegrated versions of the quark wave function and mass renormalisation factors in the on-shell renormalisation (OS) scheme. Compared to the massless case, this is a non-trivial case because the OS scheme is built at integral level and DREG leads to very simple results after integration. In any case, the unintegrated renormalisation constants are completely general and lead to a fully local cancellation of the remaining IR singularities. The treatment of UV singularities was also discussed carefully, and we constructed suitable unintegrated UV counter-terms for both self-energy and vertex corrections, which reproduce successfully the MS conditions and also lead to a fully local cancellation of UV singularities.
Finally, we tested the LTD four-dimensional unsubtraction to compute NLO QCD corrections to the decay rate of scalar and vector particles into a pair of massive quarks. The results were compared with the standard DREG expressions, and we found an impressive agreement. In particular, the transition to the massless limit is smooth because the quasi-collinear configurations of the real and virtual corrections are matched at integrand level. With the results presented in this paper, LTD four-dimensional unsubtration can be applied to any multipartonic process involving heavy quarks, and other heavy particles.  Figure 7: Total decay rate at NLO for off-shell vector particles into a pair of heavy quarks as a function of the mass, normalised to the LO. We consider three physical cases: γ * → qq (black), Z * → uū (red, up-type quarks) and Z * → dd (blue, down-type quarks). Solid lines corresponds to the results within the DREG approach, whilst the dots we obtained with the LTD unsubtracted method.

A Phase-space
The phase-space for a 1 → 2 decay with final-state particles of equal masses, p 2 i = M 2 with i = 1, 2 and s 12 the virtuality of the decaying particle, is given by with β = √ 1 − m 2 and m 2 = 4M 2 /s 12 . The corresponding 1 → 3 phase-space of the real radiation correction with an additional massless particle in the final state, (p r ) 2 = 0, is given by with where y ir = 2 p i · p r /s 12 . In order to integrate analitically the real radiation contribution, it is convenient to use the change of variables suggested in Ref. [54], i.e. y 1r = g(z) w , y 2r = g(z) z w , which allows to factorise the function h p in Eq. (101) according to In consequence, the phase-space limits, which are determined by the quadratic function h p , simplify to z ∈ [x S , x −1 S ] and w ∈ [0, 1]. The first integral in w can easily be obtained by keeping the exact -dependence. The second integral in z, however, requires to expand the result up to O( 0 ) before integration.

B Unification of coordinates
In order to avoid local mismatches in the -expansion close to the singular regions, it is necessary to unify the coordinate system and express all the on-shell momenta q i in terms of a single loop threemomentum [36]. This is a crucial point to obtain integrable expressions directly at integrand level in four space-time dimensions . Since q 3 = , we parametrise the on-shell momenta q 1 and q 2 by using the variables (ξ, v) ≡ (ξ 3 , v 3 ). Notice that each q i is set on-shell inside the corresponding dual contribution. Hence, by using the associated dispersion relations, we work with the spatial components of the momenta and fix the energy to fulfill the on-shell condition. From Eq. (16), the spatial components of the loop momenta must satisfy which leads to a system of equations, whose solution is Similarly, we can express (ξ 2 , v 2 ) in terms of (ξ, v) from q 2 = q 3 . This leads to the trivial replacement (ξ 2 , v 2 ) → (ξ, v). It is worth appreciating that this change of reference frame is well defined because the argument of the square root in Eq. (105) is always positive. In fact, due to β < 1. Then, the associated Jacobian is given by whose massless limit (i.e. β → 1) agrees with the expressions found in Refs. [36,39].

C LTD amplitudes for A * → qq(g)
In this Appendix, we collect the dual amplitudes and real squared amplitudes contributing to the NLO QCD corrections to the processes A * → qq(g), with A = {φ, γ, Z}. The tree-level vertices are given by Γ (0) φ = ı Y q 1 + c q γ 5 , Γ (0) γ = ı e e q γ µ , Γ (0) Figure 8: Momentum configuration of the NLO QCD corrections to the process A * → qq(g), under the assumption that the decaying particle does not couple to gluons.
The corresponding Born squared amplitudes, averaged over the initial-state polarisations, are The momentum configuration of the NLO QCD corrections is represented in Fig. 8.
The squared amplitudes for the real process A * → q(p 1 ) +q(p 2 ) + g(p r ), are given by where the function h p is defined in Eq. (101).