Numerical implementation of the Loop-Tree Duality method

We present a first numerical implementation of the Loop-Tree Duality (LTD) method for the direct numerical computation of multi-leg one-loop Feynman integrals. We discuss in detail the singular structure of the dual integrands and define a suitable contour deformation in the loop three-momentum space to carry out the numerical integration. Then, we apply the LTD method to the computation of ultraviolet and infrared finite integrals, and present explicit results for scalar integrals with up to five external legs (pentagons) and tensor integrals with up to six legs (hexagons). The LTD method features an excellent performance independently of the number of external legs.


Introduction
The recent discovery of the Higgs boson at the LHC represents a great success of the Standard Model (SM) of elementary particles. With the new run the primary goal is to study its properties in detail and to detect possible extentions to the SM. Precise theory predictions are needed to achieve this goal, which calls for calculations at the next-to-leading order (NLO) and beyond for multi-leg processes.
The development of automated NLO tools has seen great progress in recent years. Computing higher-order corrections in QFT, in particular in QCD and in the EW sector of the SM is highly challenging. The complexity increases as the number of external particles gets bigger and the order of the perturbative expansion. The task is far from trivial and each step presents its own difficulties: one needs first to generate the virtual and real scattering amplitudes, then carry out the integration over the loop momenta -1 -for the virtual contribution and finally perform the phase-space integration for both real and virtual corrections after taking proper care so that the infrared divergencies cancel. In particular, infrared singularities of the virtual contribution can be subtracted by using appropriate semi-analytical terms and combine them with the ones stemming from the real corrections to produce finite results [1]. Purely numerical approaches on the integration of loop momenta have been discussed extensively in the literature [2][3][4][5][6][7][8][9][10][11][12][13][14]. The generation of amplitudes and calculation of cross-sections at one loop has seen great progress in recent years and algorithmic calculations at NLO are now considered standardised, based on purely numerical [15][16][17] and a mix of analytical and numerical approaches [18][19][20]. Substantial progress has also been made at higher orders [21][22][23].
The loop-tree duality (LTD) method [24][25][26][27][28][29][30][31][32][33] establishes that generic loop quantities (loop integrals and scattering amplitudes) in any relativistic, local and unitary field theory can be written as a sum of tree-level-like objects obtained after making all possible cuts to the internal lines of the corresponding Feynman diagrams, with one single cut per loop and integrated over a measure that closely resembles the phasespace of the corresponding real corrections [24,25]. This duality relation is realised by a modification of the customary +i0 prescription of the Feynman propagators and encodes the causal structure of the scattering amplitudes in the expected way. The analysis of the singular behaviour of one-loop integrals and scattering amplitudes in this framework at the integrand level in the loop momentum space shows that there is a partial cancellation of singularities among different dual contributions such that physical infrared and threshold singularities remain restricted to a compact region of the loop three-momentum [28][29][30]. This feature opens up the intriguing possibility that virtual and real radiative corrections can be brought together under a common integral and be treated simultaneously with Monte Carlo techniques though a convenient mapping of the external momenta entering the virtual and real scattering amplitudes [33].
In this work, we present a first numerical implementation of the LTD method and we apply it to the computation of multi-leg one-loop scalar and tensor integrals. The outline of the paper is as follows. In Section 2 we review the LTD method at one-loop and discuss the singular behaviour of the dual integrand in the loop momentum space. In Section 3 we introduce the contour deformation in the loop three-momentum space which is used for the numerical loop integration. We present explicit numerical results for various external momenta configurations, in Section 4 for scalar integrals up to pentagons and in Section 5 for up to rank 3 tensor integrals with five and six external legs. Finally, our conclusions are presented in Section 6.
- 2 -The FTT and the duality relation can be illustrated with no loss of generality by considering their application to the basic ingredient of any one-loop Feynman diagrams, namely a generic one-loop scalar integral L (N ) with N (N ≥ 2) external legs. The momenta of the external legs are denoted by p µ 1 , p µ 2 , . . . , p µ N and are clockwise ordered ( Fig. 1). All are taken as outgoing. To simplify the notation and the presentation, we also limit ourselves in the beginning to considering massless internal lines only. Thus, the one-loop integral L (N ) can in general be expressed as: where q µ is the loop momentum (which flows anti-clockwise). The momenta of the internal lines are denoted by q µ i ; they are given by and momentum conservation results in the constraint The value of the label i of the external momenta is defined modulo N, i.e. p N +i ≡ p i .
The number of space-time dimensions is denoted by d (the convention for the Lorentzindices adopted here is µ = 0, 1, . . . , d − 1) with metric tensor g µν = diag(+1, −1, . . . , −1). The space-time coordinates of any momentum k µ are denoted as k µ = (k 0 , k), where k 0 is the energy (time component) of k µ . It is also convenient to introduce light-cone coordinates Throughout the paper we consider loop integrals and phase-space integrals. If the integrals are ultraviolet or infrared divergent, we always assume that they are regularized by using analytic continuation in the number of space-time dimensions (dimensional regularization). Therefore, d is not fixed and does not necessarily have integer value.

Loop-tree duality at one-loop
We consider a general one-loop N -leg scalar integral (see Fig. 1) in dimensional regularization, with d the number of space-time dimensions, are Feynman propagators that depend on the loop momentum , which flows anticlockwise, and the four-momenta of the external legs p i , i ∈ α 1 = {1, 2, . . . N }, which are taken as outgoing and are clockwise ordered. The momenta of the internal lines are denoted as q i,µ = (q i,0 , q i ), where q i,0 is the energy (time component) and q i are the spatial components, which are defined as q i = + k i with k i = p 1 + . . . + p i , and k N = 0 by momentum conservation. We also define k ji = q j − q i which, in fact, is independent of the loop momentum . The corresponding dual representation of the scalar integral in Eq. (2.1) is obtained from the loop-tree duality (LTD) theorem [24]: -3 -are dual propagators, η is an arbitrary future-like vector, i.e., a d-dimensional vector that can be either light-like (η 2 = 0) or time-like (η 2 > 0) with positive definite energy η 0 ≥ 0, and selects the internal loop on-shell modes, G −1 F (q i ) = 0, with positive definite energy, q i,0 ≥ 0. Hence, the LTD theorem expresses the usual loop Feynman integral, Eq. (2.1), as a sum of single-cut phase-space integrals, Eq. (2.3), with (2.6) as the single-particle phase-space integration measure. The LTD theorem is valid not only for scalar one-loop integrals, but can straightforward be extended to deal with scattering amplitudes [24] and higher orders of the perturbative expansion [26,27]. The integrand of the dual representation of one-loop integrals or scattering amplitudes feature certain types of singularities leading to ultraviolet, infrared or threshold singularities. This singular behaviour has already been thoroughly discussed in [29,33]. We briefly recapitulate here the main points that are relevant in the present context. are the so-called dual propagators, as defined in Ref. [20], with η a future-like vector, η 2 ≥ 0, with positive definite energy η 0 > 0. The delta functionδ (q i ) ≡ 2π i θ(q i,0 ) δ(q 2 i − m 2 i ) sets the internal lines on-shell by selecting the pole of the propagators with positive energy q i,0 and negative imaginary part. In the following we take η µ = (1, 0), and thus −i0 η k ji = −i0 k ji,0 . This is equivalent to performing the loop integration along the on-shell forward hyperboloids. Let us mention that in the light-cone coordinates (ℓ + , ℓ − , l ⊥ ), where ℓ ± = (ℓ 0 ± ℓ d−1 )/ √ 2, Feynman propagators vanish at hyperboloids in the plane (ℓ + ,ℓ − ) which are similar to those depicted in Fig. 1 but rotated by 45 degrees. Consequently, by selecting the forward hyperboloids the integration limits of either ℓ + or ℓ − are restricted and the restrictions are different for each dual integral. For this reason, although Eq. (3) is valid for any system of coordinates, we will stick for the rest of the paper to Cartesian coordinates where all the dual integrals Figure 2: On-shell hyperboloids for three arbitrary propagators in Cartesian coordinates in the ( 0 , z ) space (left). Kinematical configuration with massless propagators leading to infrared singularities (right). In the latter case, the on-shell hyperboloids degenerate to light-cones.
For generic masses, the loop integrand in Eq. (2.1) becomes singular at the on-shell hyperboloids defined by q . This is illustrated in Fig. 2 for a given kinematical configuration with three internal loop propagators. Solid lines in Fig. 2 represent the forward on-shell hyperboloids, and dashed lines the backward on-shell hyperboloids. The LTD method is equivalent to evaluating the sum of the integrals along the forward on-shell hyperboloids with the singularities appearing at the intersection of each forward on-shell hyperboloid with the forward of backward on-shell hyperboloid of the other propagators. A crucial point of this discussion is the observation that dual propagators can be rewritten as i,0 can be interpreted as the loop energy measured along the forward on-shell hyperboloid with origin at −k i . From Eq. (2.7) it is obvious that dual propagators become singular, G −1 D (q i ; q j ) = 0, if one of the following conditions is fulfilled: The first condition, Eq. (2.8), is satisfied if the forward on-shell hyperboloid of G F (q i ) intersects with the backward on-shell hyperboloid of G F (q j ). The second condition, Eq. (2.9), is true when the two forward on-shell hyperboloids intersect each other. The solution to Eq. (2.8) is an ellipsoid in the loop three-momentum space and requires k ji,0 < 0. Moreover, since it is the result of the intersection of a forward with a backward on-shell hyperboloid the distance between the two propagators has to be future-like, k 2 ji ≥ 0. Actually, internal masses restrict this condition to The second equation, Eq. (2.9), leads to a hyperboloid in the loop three-momentum space, and there are solutions for k ji,0 either positive or negative, namely when either of the two momenta is set on-shell. Here, the distance between the momenta of the propagators has to be space-like, although also time-like configurations can fulfill Eq. (2.9) as far as the time-like distance is small or close to light-like: As it was demonstrated in [29], the integrand singularities appearing from the intersection of forward with forward on-shell hyperboloids cancel among dual contributions. To see that one needs to keep in mind that propagators are positive inside the on-shell -5 -hyperboloids and negative outside. When integrating along the forward on-shell hyperboloids, every singularity is crossed twice. Firstly when going from the inside to the outside (or from the outside to the inside) and secondly from the outside to the inside (or from the inside to the outside). The crucial point is that the contributions coming from the two integrands have opposite sign and thus cancel out. Note that the imaginary dual prescription η ·k ji changes sign from the one dual contribution to the other to ensure the cancellation of the singularities. On the contrary, the singularities from the intersection of a forward with a backward on-shell hyperboloid survive because only a single dual contribution leads to that singularity and there is no possibility of cancellation. In the case of integrable singularities, a contour deformation can be employed as explained in the next section.
The action of the LTD can be encoded symbolically by the following matrix scheme Each line in the matrix to the right of the arrow in Eq. (2.12) represents a dual contribution with one single propagator on-shell, δ = δ(q i ). The column index points to the corresponding dual propagators, G D = G D (q i ; q j ). This scheme can now be used to graphically indicate the position of different singularities in a given dual integral. In Figs. 3 and 4 we apply it to a triangle and a box respectively. To be more specific, in Fig. 3 each of the 3D plots in the r.h.s represents the singularities of any one dual contribution. We plot the ellipsoid (orange surfaces) and hyperboloid (blue surfaces) singularities in the loop three-momentum space. The blue dots are the foci of the onshell hyperboloids, i.e. −k i , i ∈ {1, 2, 3}. In the l.h.s, we see the singularity scheme, where the first line of the matrix corresponds to the first plot in the r.h.s, the second line corresponds to the second 3D plot and so on. In the matrix, an H indicates that the corresponding dual propagator from Eq. (2.12) generates an hyperboloid singularity, E stands for ellipsoid singularities, and zero means no singularity. Similarly, for a four-point function in Fig. 4. In both cases, the hyperboloid singularities always appear pairwise across the dual contributions. This is not by accident. Due to the symmetry of Eq. (2.9) under the exchange of i (i counts dual contributions) and j (j counts leg positions) the hyperboloid singularities always appear in pairs and are distributed symmetrically around the main diagonal. Inspecting Eq. (2.8), which is the defining equation for ellipsoid singularities we see that this equation is not symmetric under the exchange of indices. Thus for every ellipsoid singularity we have a zero as its counterpart.  At this point, we have established that the hyperboloid singularities do cancel among dual contributions and therefore we do not need to treat them in any special manner. Still though, they do have an impact on the way we need to deform our contour. This is due to the fact that in order to preserve the cancellation of hyperboloid singularities, dual contributions featuring the same hyperboloid singularity must receive the same deformation. To further illustrate this point, let us look at the pentagon example shown in Fig. 5. In Fig. 5, contributions one, two and three are coupled via their common hyperboloid singularities. Thus, they need to receive the very same deformation that accounts for all the ellipsoid singularities occurring within those contributions. These are found at position four of the second contribution and positions one and four of the third contribution. The fourth dual contribution is not coupled to any other contribution and a standalone deformation can be applied. The fifth contribution does not require any treatment.
Contributions are coupled: Every contribution receives all deformations that occur within the group. → Deform with ellipsoids that itself contains. → No deformation needed here. As a general strategy, one organises the dual contributions into groups. A group is a set of pairwise coupled contributions. To each of the groups a contour deformation is applied independently from the others. Within a group every contribution receives the same deformation that accounts for all the ellipsoids of the group. Turning back to the example in Fig. 5, we have three groups: the first group involves contributions one to three, the second group is contribution four and the third group is contribution five.

The deformation of the contour
As we saw in Section 2, the ellipsoid singularities (forward-backward type) lead to integrable threshold singularities that lie on the real axis. To deal with them, we need to deform the integration path into the imaginary space. Every valid deformation must satisfy a certain set of requirements [3]: 1. The deformation has to respect the i0-prescription of the propagator. In general, a contour deformation in the loop three-momentum space has the form: where κ is a function of the loop momentum and the external momenta. In our case, we want to perform the integration over a product of dual propagators. Inserting Eq. (3.1) into the on-shell energy relation, we obtain The i0-prescription tells us in which direction to deform when coming close to a singularity. Hence, any valid deformation must match this prescription. Consequently we need to have The deformation should vanish at infinity: We are looking for a deformation that does not change the actual value of the integral. We do not want |κ| to grow for | | → ∞. An easy way to satisfy this condition is to choose κ such that |κ| → 0 as | | → ∞. 1 With these conditions in mind, we construct the deformation in the following way: As explained in Section 2, we first organise the dual contributions into groups. For every ellipsoid singularity of the group we include a factor: with q i = + k i and the loop three-momentum. The deformation factor in Eq. (3.4) consists of two main components. The first component defines the direction of the deformation, and is given by the sum of the two unit vectors q i / q 2 i and q j / q 2 j . As shown in Fig. 6 the vectors q i and q j have their origin in −k i and −k j , respectively, and the deformation is designed to point to the outside of the singularity ellipsoid. For an efficient numerical implementation, however, we should also take into account in the selection of the deformation parameters that for massive propagators the vectors −k i and −k j might be slightly displaced from the true focal points of the ellipsoid. Inside the ellipsoid, the sum of the two unit vectors q i / q 2 i and q j / q 2 j helps to flatten the deformation and indeed they cancel each other along the major axis of the ellipsoid. By choosing all the scaling parameters λ ij < 0 for all possible combinations {ij} we satisfy the first condition in Eq. (3.3).
The second component in Eq. (3.4) is the exponential factor exp(−G −2 D (q i ; q j )/A ij ) that suppresses the deformation at infinity. At singular points, G −2 D (q j ; q i ) vanishes and thus the deformation reaches its maximum. Far away form the singularity, −G −2 D (q j ; q i ) reaches a large negative value and thus the exponential tends rapidly to zero. Finally, the factor λ ij is a scaling factor, and A ij determines the width of the deformation. The indices ij in λ ij and A ij indicate that those parameters can be selected individually for each deformation contribution for optimization purposes. Then, we sum over the entire group of coupled singularities and arrive at: The corresponding Jacobian can be calculated analytically or numerically; in our current implementation we have chosen the analytic way.

Results for multi-leg scalar one-loop integrals
We have implemented the LTD method in a C++ code and all the results in this paper were obtained on a desktop machine with an Intel i7 (3.4GHz) processor with 8 cores and 16 GB of RAM. The program uses the Cuba library [35] as a numerical integrator. The user needs only to input the number of external legs, the external momenta, the internal masses and has the freedom to change the parameters of the contour deformation. The momenta and masses can be read in from a text file. The user can choose between Cuhre [36,37] and VEGAS [38], and give the desired number of evaluations. At run time, the code performs the following steps: 1. Reads in masses and external momenta.
3. Checks where hyperboloid singularities occur, groups the dual contributions accordingly and applies the contour deformation.

Calls the integrator.
We use MATHEMATICA 10.0 [39] to generate random momenta and masses to scan as much of the phase-space as possible, to ensure that the program works properly in all regions. For our numerical results, the routine Cuhre was used unless otherwise stated. The momenta and masses of all the sample phase-space points used in the following sections are collected in Appendix A. We mainly used LoopTools 2.10 [34] and also SecDec 3.0 [40] to produce reference results to compare with.

Scalar triangles
We consider first infrared finite scalar triangle integrals. The sample point P1 in Table 1 has all internal masses equal while P2 has three different internal masses. Momenta and masses were chosen randomly between −100 GeV and +100 GeV. Similarly, P3 in Table 1 has all internal masses equal whereas in P4 all three of them have different values.
Points with momentum configurations that do not need deformation (i.e. whose loop integral is purely real) are computed in well below one second with a precision of at least 4 digits. For points with momentum configurations that require deformation, the calculation time increases to typically 3 − 15 seconds.

Scalar Triangle Real Part
Imaginary Part P1 LoopTools  An important check of our implementation is the mass-scan around threshold. In Fig. 7

Scalar Boxes
Next we consider infrared finite box scalar integrals. To get good precision (4 digits), for boxes that need deformation we use 4 to 5 · 10 6 Cuhre calls, whereas for phase-space points with no deformation only 5·10 4 calls, the same as in the triangle case. This is reflected in the running times. Points with deformation require about 15 seconds whereas points with no deformation well below one second. While it is practically guaranteed to get the no-deformation points with good precision, for points with deformation the quality of results depends on the proper choice of the deformation parameters. Therefore, we mainly focus our attention to such points in the following. The sample points P5 and P7 of Table 2 correspond to a momentum configuration in which all four internal masses are equal. In P6 and P8 all masses are different. In P9 two adjacent internal lines have equal masses as well as the two opposing ones. P10 represents a situation in which opposite lines have equal masses. We perform again a mass-scan (see Fig. 8) with all internal masses equal, i.e. m i = m, i ∈ 1, 2, 3, 4. The center-of-mass energy s was kept constant while the mass m was varied. The program deals well with all kinds of boxes, even when many different kinematical scales are involved. In Fig. 8, two thresholds are crossed at 2m/ √ s = 0.65 and 1. From right to left, the number of ellipsoid singularities grows by one after each threshold is crossed, starting from one to end up to three.

Scalar Pentagons
Let us now turn to pentagon diagrams. No-deformation points are computed with 10 5 evaluations which takes about 0.5 seconds. Points with deformation demand 5 · 10 6 evaluations to maintain the level of precision of the triangles and boxes. This results to an average calculation time of about 30 seconds. In Table 3 we display a collection of pentagon example results for different kinematical configurations. In P11 and P13 all internal masses are equal; in P14 they Our implementation of the LTD method shows its robustness by producing accurate results regardless of the kinematical situation. This statement is further supported by an energy-scan which we performed and which is shown in Fig. 9. The center-of-mass energy s is varied. This is realized by varying p 3 while keeping p 2 3 constant. Of course, due to momentum conservation, this involves p 2 4 = (p 1 + p 2 + p 3 ) 2 not being constant. In this scan, we cross three thresholds at s ≈ −8.5 · 10 3 , −13.5 · 10 3 and −21 · 10 3 GeV 2 which divide the scan into four zones. From right to left, we start with zero ellipsoid singularities in the first zone, then we have one in the second zone, two in the third zone and finally one in the leftmost zone.

Tensor loop integrals
The LTD relation for scalar loop integrals can easily be extended to deal tensor integrals. As long as the quantum field theory is local and unitary, these tensor factors do not lead to additional singularities [24] and the LTD method can then be applied in a straightforward manner. If the one-loop integral features a non-trivial numerator N ( , {p i }), then, the LTD theorem takes the form While the numerator is formally left unchanged, there is actually a potential implication. The presence of the dual delta function demands q In other words, whenever we perform a single cut of a Feynman graph, the numerator has to be evaluated at the position of the cut which is fixed by the dual delta function. As a direct consequence, the numerator takes a different form in each dual contribution. Another important aspect to take into consideration is the cancellation of singularities among dual contributions. Here, we would like to make explicit that the numerators do not spoil the cancellation of the hyperboloid singularities. A typical numerator is a polynomial of scalar products of the loop-momentum with external momenta: · p k .
Let us see what happens to a single factor when it hits the singularity. Note first, that the hyperboloid singularity is given by Eq. (2.9) which we rewrite in the more suitable form Using Eq. (5.3), the loop-momentum contracted with some external momentum p k is: where we have used Eq. (5.4). This means that the numerators of two dual contributions i and j take the same value at their common pole, thus leaving the cancellation of -15 -hyperboloid singularities intact. This is an important property to take advantage of, because it allows us to straightforwardly apply the LTD method to such diagrams without any additional effort.

Tensor Pentagons
Next, we investigate tensor pentagon integrals at the one-loop level with numerators up to rank three. The number of evaluations is chosen to be the same as in the scalar case, i.e. 10 5 times for no-deformation points and 5 · 10 6 times for phase-space points that require deformation. This results in calculation times of about 1 second and 30 seconds, respectively. In Table 4 we show a selection of sample points. The reference points P16 and P18 feature the rank two numerator ( · p 3 ) × ( · p 4 ) while P17 and P19 have the numerator ( · p 3 ) × ( · p 4 ) × ( · p 5 ). All the points have all internal masses equal. P19 actually contains six ellipsoid singularities whereas the other sample points have two to three. We include this point to demonstrate that the program does well even under such challenging circumstances.
For tensor pentagons and hexagons, we have used the program SecDec [40] to crosscheck our results. We have run SecDec taking no care to optimise its runtime. This means that in the following, whenever we present the running times of SecDec we do it for completeness reasons and not because we imply that our code compares better or worse to SecDec. A proper comparison of our implementation with available codes is beyond the scope of this paper.
We have performed several different scans; a sample is presented in Fig. 10. In this energy-scan, we varied p 1 and thus the center-of-mass energy s = (p 1 + p 2 ) 2 , similar to what we have done with scalar pentagons. The corresponding numerator function is ( · p 1 ) × ( · p 2 ) × ( · p 3 ), which means that both numerator and denominator change in the scan. In Fig. 10, one can see that the LTD method is able to successfully pass this test.

Tensor Hexagons
Finally, we compute hexagon tensor integrals. The number of evaluations for nodeformation points is 10 6 and for deformation points 8 · 10 6 . The typical corresponding computation times are about 10 and 75 seconds, respectively. Since LoopTools can provide results only up to pentagons, we used exclusively here the program SecDec to cross-check.
We present a selection of sample points in Table 5. P20 and P22 feature the rank one numerator · p 1 , the former has all internal masses different, in the latter they are all equal. P21 has six distinct internal masses and the numerator function   Figure 10: Energy-scan of a rank three tensor pentagon around the threshold region. The red curve is done with LoopTools and the blue points are obtained with the LTD method.

Conclusions
The Loop-Tree Duality has many appealing theoretical properties for the calculation of processes with many external legs. In this paper, we have investigated the practicability of a first numerical implementation of the LTD method. In our analysis of the singular behaviour of the loop integrand, we found two distinct types of singularities: Ellipsoid singularities which require the application of contour deformation and hyperboloid singularities that occur pairwise and cancel among different dual contributions. In order to preserve their cancellation, dual contributions featuring the same hyperboloid singularity pair must receive the same contour deformation. This leads to the following algorithm: Sets of pairwise coupled dual contributions are organized into groups. Each group is deformed independently from the others and each dual contribution of such a group receives the exact same contour deformation that accounts for all the ellipsoid singularities of the entire group.
We applied a contour deformation that efficiently deals with the ellipsoid singularities by meeting all the important criteria [3]. This setup has been successful in the calculation of finite multi-leg scalar and tensor integrals. We found the results to be in very good agreement with the reference values produced by LoopTools and SecDec. An important further check of our implementation presented here was various scans which show that the code handles equally well broad slices of the phase space. The code excels in cases that involve many external legs as it shows a modest increase in running times in comparison to cases with fewer legs. From this first study, we can be optimistic that our implementation of the LTD method offers a competitive alternative -18 -for computing multi-scale, multi-leg scalar and tensor one-loop integrals.
In this paper we have only considered IR-and UV-finite integrals. However, the treatment of IR-and UV-divergent graphs in the context of the LTD method and the description of how to combine directly the virtual with the real radiative corrections in order to obtain an infrared finite implementation has been presented recently [33]. The extension of our code to deal with these cases is an ongoing project.