A stroll through the loop-tree duality

The Loop-Tree Duality (LTD) theorem is an innovative technique to deal with multi-loop scattering amplitudes, leading to integrand-level representations over an Euclidean space. In this article, we review the last developments concerning this framework, focusing on the manifestly causal representation of multi-loop Feynman integrals and scattering amplitudes, and the definition of dual local counter-terms to cancel infrared singularities.


Introduction
The most successful description of the microscopic structure of Nature is currently given by the Standard Model (SM), which is based on quantum field theories (QFT) with specific gauge symmetries. The strong interplay between experiment and theory relies on confronting precise predictions with accurate data. In the context of high-energy particle physics, the Large Hadron Collider (LHC) and the planned future colliders [1][2][3][4][5][6][7][8] will keep on reducing the experimental uncertainties, forcing to reach a compatible accuracy level from the theory side. However, there are huge bottlenecks that prevent straightforward calculations.
Due to their complexity, exact solutions of QFTs are not available for generic scattering processes. Thus, extracting predictions from theory requires to use approximations, whose applicability is restricted to specific situations. Whilst lattice methods are reliable in the low-energy regime, perturbation theory is the most suitable strategy to tackle the phenomenological description of particle collisions at high-energies. Even if the perturbative approach reduces the original calculation to a series expansion in the interaction couplings, adding higher-orders is far from trivial due to the presence of complicated phase-space and Feynman loop integrals. Moreover, these objects feature numerical instabilities originated in threshold configurations, as well as other singularities that cancel only after putting all the contributions together.
In the last forty years, a huge effort has been done to develop new methodologies for an efficient calculation of physical observables at higher-orders. As a first step, proper regularization methods have been applied in order to make explicit the singular structures of these objects. One of the customary choices is Dimensional Regularization (DREG), which consists in modifying the number of space-time dimensions to achieve integrability [9][10][11][12]. However, in the context of QFT, changing the number of space-time dimensions leads to definition problems, such as the γ 5 [13,14], and prevents a straightforward numerical implementation. For this and other technical reasons, there is an ongoing effort in the highenergy physics community to develop new strategies that locally regularise QFT whilst keeping the standard four dimensions of the space-time [15][16][17][18][19][20].
On the other hand, loop contributions must be combined with real-emission corrections and suitable counter-terms to cancel the infrared (IR) and ultraviolet (UV) singularities, leading to a reliable prediction of the physical cross-sections and differential distributions. There are several strategies to achieve this cancellation, most of them based on the subtraction methods [55][56][57]. Generally speaking, these frameworks involve adding and subtracting counter-terms which locally cancel the IR behaviour of the real part, but at the same time they are easy enough to be integrated analytically and subtracted from the virtual contribution. There are several variations of this approach: dipole subtraction [57,58], q T subtraction/resummation [59,60], antenna subtraction [61], colorful subtraction [62,63], among others. Also, there are alternatives that rely on a different way to achieve the cancellation of singularities, like the n-jetiness [64,65] or the local analytic subtraction frameworks [66,67].
Besides those methods, a novel strategy to tackle, simultaneously, the efficient calculation of loop integrals and physical observables was developed. This framework is based on the Loop-Tree Duality (LTD) theorem [68][69][70][71], which opens loops into trees and recasts virtual states into configurations that resemble real-radiation processes. From the mathematical point of view, the LTD transforms the integration domain of loop integrals into an Euclidean space. In this way, a more intuitive understanding of the regions responsible for the singular structure of the loop integrals emerges [72]. This knowledge can be used to explore novel techniques pointing towards more efficient numerical implementations [73,74], integrand-level simplifications through asymptotic expansions [75][76][77] and unveiling the universal structures of scattering amplitudes at higher-orders [78][79][80][81]. Moreover, the LTD provides the perfect framework to handle cross-section calculations in a fully unified way. Since the dual representation of loop integrals is defined in Euclidean domains, the virtual and real contributions can be directly combined at integrand-level to achieve a fully local cancellation of IR singularities. This is the so-called Four-Dimensional Unsubtraction (FDU) approach [82][83][84][85], that profits from a local cancellation of singularities which makes it possible to bypass additional regulators (such as DREG). Also, this formalism allow us to write local UV counter-terms that exactly reproduce the expected results in conventional renormalisation schemes [80,81,84,85], as well as fully local IR dual counter-terms [86].
During the last two years, we have explored new features of the LTD approach. In particular, we discovered that a more general analysis of the singular structure of multi-loop multi-leg amplitudes is possible. By inquiring into the physical and anomalous thresholds of one and two-loop amplitudes [87], we proved that it was possible to get rid of non-physical singularities and retain only those contributions compatible with causality. Several studies about the causal structure of scattering amplitudes are available in the literature, using different techniques [88][89][90][91]. In Ref. [92], we presented for the first time a manifestly causal integrand-level representation inspired by the LTD theorem. This was applied to remove unphysical threshold singularities and obtain causal integrand-level representations of several topological families of multi-loop multi-leg Feynman integrals [92][93][94][95][96][97][98][99]. Moreover, all-order causal formulae were obtained using novel algebraic relations [100], leading to an automatised framework to implement these calculations [101]. A purely geometrical interpretation [102] was also developed, showing a complementary approach to tackle the causal structure of multi-loop multi-leg scattering amplitudes and Feynman integrals.
The purpose of this review is to give a brief summary of the latest developments regarding the LTD-based methods. To ease the presentation, it is organised in two parts. In the first one, starting in Sec. 2, we settle the notation of the LTD formalism. Special emphasis is put on the manifestly causal representation of loop integrands in the Euclidean space, as well as its derivation from Cauchy's residue theorem. The basis for a rigorous mathematical treatment of the nested residue strategy is given in Sec. 3. Then, we discuss the advantages of the Euclidean representation of loop integrands. In Sec. 4, we focus on the asymptotic expansion, whilst we discuss how to obtain compact causal formulae for generic topological families of multi-loop multi-leg Feynman integrals and scattering amplitudes in Secs. 5 and 6. In Sec. 7, a few words are given to explain novel developments suggesting a deep connection among graph theory, algebra and the causal structure of multi-loop scattering amplitudes.
After that, we start the second part of the review, where we center into the computation of physical observables exploiting the LTD approach. In Sec. 8 we discuss the basis of the FDU formalism, showing a physical examples. Then, in Sec. 9, the construction of local counter-terms for cancelling IR singularities is presented with more general examples, including γ * → 3 jets at NLO. Finally, we present the conclusions and future research directions in Sec. 10.

Causality within the LTD formulation
Among the several advantages of the LTD formalism, we have recently explored the manifestly causal representations. This constitutes a big step towards an efficient numerical implementation of this approach, since spurious non-physical and non-causal singularities are completely avoided. In the following, we will introduce some useful notation and describe the basis of the formalism. In particular, we will explain how the LTD computation through the nested residue strategy leads to the causal representations of different topological families to all-loop orders.

Dual scattering amplitudes
A generic L-loop scattering amplitude with N external legs {p j } N , and n sets of internal lines, each set being defined by the specific dependence on the loop momenta, is written as This corresponds to an integral in the Minkoswki space of the L-loop momenta, { s } L , involving the product of Feynman propagators, G F (q i ) = (q 2 i − m 2 i + ı0) −1 , and numerators, N ({ s } L , {p j } N ), given by the Feynman rules of the theory considered. The d-dimensional integration measure in DREG reads and the usual Feynman propagator for a single particle is is the positive on-shell energy of the loop momentum q is , written in terms of its spatial components q is , its mass m is , and the Feynman's infinitesimal imaginary prescription ı0. Then, we introduce the shorthand notation to express the product of Feynman propagators of one set or the union of several sets. In the former formula, s represents the set of all internal propagators with internal momenta q is = s + k is , that depend on the loop momentum, s , or a specific linear combination of loop momenta, and a combination of external momenta k is , with i s ∈ s and a is arbitrary powers. The LTD representation for scattering amplitudes is obtained by the iterative application of the Cauchy's residue theorem (CRT) integrating out one degree of freedom for each loop momentum, and closing the Cauchy contour from below the real axis, selecting the poles with negative imaginary part in the complex plane of the loop momentum. This results in the modification of the infinitesimal complex prescription of the Feynman propagators [68,69], that needs to be considered carefully in order to preserve causality. Then, starting from Eq. (2.1), we set on-shell the propagators that depend on the first loop momentum, q i 1 , and define where taking the residue is equivalent to integrating out the energy component of the loop momenta. Now, we construct the nested residue iterating until the r-th set as, All sets before the semicolon contain one propagator that has been set on shell and are linearly independent, while all the remaining propagators are kept off-shell. Thus, this representation is equivalent to open the loop amplitude to non-disjoint trees. Finally, the integration measure after integrating out the energy component is modified according to transforming the d-dimensional Minkowski space into a (d − 1)-dimensional Euclidean one.

Multi-loop topologies through the LTD
By the properties obtained from the nested residues, Eq. (2.8), we are able to construct multi-loop amplitudes from the LTD representation which are compact and manifestly causal to all orders. The first topology under consideration is called Maximal Loop Topology (MLT), and is characterised by L + 1 momentum sets where the momenta of the first L sets depend on one single loop momentum The momenta of the extra set, L + 1, are given by a linear combination of all the loop momenta, namely q i L+1 = − L s=1 s + k i L+1 . Here, k is and k i L+1 denote linear combinations of external momenta. The LTD representation of the MLT topology, displayed in Fig. 1a, is presented in a compact and symmetric form by evaluating the nested residue of Eq. (2.8), which leads to where the bars within the integrand indicate a reversal of momentum flow, qī s = −q is , which is needed to preserve causality, and is equivalent to select the on-shell modes with negative energy of the original momentum flow. Each term of the sum in the integrand of Eq. (2.11) contains one set i with all its propagators off-shell, while the remaining L sets contain a single on-shell propagator each; a necessary condition to open multi-loop amplitudes into disjoint trees. The dual representation of Eq. (2.11) becomes singular when one or more off-shell propagators eventually become on-shell, generating a disjoint tree dual sub-amplitude. It is worth mentioning that the LTD representation exhibits an interesting structure when all the contributions are added together. For example, in the case of a MLT configuration with one propagator in each loop set and one incoming and outgoing momentum, we obtain A (L) This expression is free of unphysical singularities, and written in terms of the so-called causal propagators, λ ± 1 . The causal propagators encode all the possible physical singularities that might occur. Both λ ± 1 in the integrand of Eq. (2.12) are associated to physical thresholds, as shown in Fig. 2. Additionally, Eq. (2.12) is independent of the initial momentum configuration in the Feynman representation. Going a step further in the topology complexity, we consider a topology containing one extra set of momenta that depends on the sum of two loop momenta, q i 12 = −( 1 + 2 )+k i 12 , denoted as 12. This configuration, as can be appreciated in Fig. 1b, is called Next-to-Maximal Loop Topology (NMLT), characterised by L + 2 sets of propagators, with each set categorised by the dependence on a specific loop momentum or a linear combination of the L independent loop momenta. The LTD representation for the NMLT is given by the compact and factorised expression     where the sets after the semicolon are put off shell. Following the previous procedure, we also presented in Ref. [92] the LTD representation for the Next-to-Next-to-Maximal Loop Topology (N 2 MLT), characterised by L + 3 sets of propagators and depicted in Fig. 1c. Again, factorised formulae for the dual representation are obtained, namely    where the last two terms of the rhs in the integrand remain fixed by the condition that the sets (2,3,23) cannot generate a disjoint subtree. We can observe that the second term in Eq. (2.17), contain a two-loop subtopology involving five sets of momenta grouped into three sets. Therefore, the propagators in the set 1 and 23 cannot be simultaneously off shell. Noteworthy, there are very compact explicit formulae for the NMLT and N 2 MLT configurations which make use of the causal propagators. In Sec. 5, we will enter into more details, and provide a nice conceptual interpretation of the manifestly causal dual representation in terms of entangled thresholds.

Mathematical properties of the nested residues
Generic scattering amplitudes are defined by integrals of rational functions. As already explained in Sec. 2, the multi-loop LTD framework is based on the CRT and, in this section, the formal foundations for an L-loop Feynman diagram and some of their immediate consequences are explained. As it was discussed in Ref. [97], iterated residues can be easily computed taking into account the quadratic structure of Feynman propagators. Let us start with the space C (R L ) of the functions with domain R L and co-domain C. If the natural inclusion of R into C is denoted by i, then the iterated residues are defined as the recursive application of the functor −Res • i : C (R n ) → C (R n−1 ) , where Res represents the residue of the argument for all its negative imaginary part poles and the minus sign appears in agreement with CRT as the integration is performed clockwise.
For the computation of multi-loop Feynman integrals and scattering amplitudes, internal momenta are given as linear combinations of external and loop momenta. For instance, for the scalar sunrise diagram, the associated integrand has the form In this case, the first iteration of the iterated residues can be developed with respect to the variable q 1,0 . This means that we extend q 1,0 to the complex plane, and as q 2,0 is still a real variable, the poles are located in the complex plane as shown in Fig. 3, where the dashed line represents the poles ±q (+) 3,0 − q 2,0 − k 3,0 . As they depend on q 2,0 , they can be located at some point along each of the dashed lines. Then, after the computation of the first iterated residues, we obtain The function obtained after the first residue iteration is given as a sum of terms. The next step corresponds to the extension of the variable q 2,0 to the complex plane. There, the pole structure of each term of Eq. (3.2) is quite similar. This can be represented pictorially as shown in the Fig. 4, where the pole structure of each term is presented as a complex plane term. The gray blob in both terms represents the dependence of the imaginary part of the pole q (+) 2,0 − k 3,0 on the three-momenta q 2 and q 3 , having no definite imaginary part sign. These poles having no definite imaginary part sign are, in general, called displaced poles, and they happen to cancel their contributions to the iterated residues through the relation for a i and a ij linear combinations of energy components of external momenta. Due to this cancellation, displaced poles can be ignored in the computation, leading to the concept of nested residues. Furthermore, the final result of the nested residues is independent of the order of the iteration, although the expressions are not identical term by term. A rigorous proof of Eq. (3.3) is shown in Ref. [97]. As a consequence of the nested residue strategy, we can infer a general formula for the causal structure of the MLT diagram (see Fig. 1a). This is directly reached applying partial fractions to each term in the nested residue, leading to Additionally, this result can be used for the insertion of an MLT diagram within a more general topology. In this manner, it is possible to consider any MLT insertion as a single propagator whose on-shell energy is the sum of the on-shell energies of all the internal propagators of the MLT. Furthermore, with the computation of the nested residues, every family of a given topological complexity k, N k−1 MLT, can be re-expressed as a linear combination of convolutions of lower topological complexity families, as it is for the NMLT and N 2 MLT, both defined by their amplitudes, This decomposition is shown in Figs. 5 and 6 for each of these topological families. It is worth mentioning that the convolution symbol does not represent a pure factorisation, as it implies the use of the on-shell conditions to express all the off-shell momenta. The NMLT(L) diagram decomposes into two terms with two MLT diagrams. One of these terms contains a MLT(L − 2) in convolution with a MLT(2), while the other term contains a convolution of a MLT(L − 2) diagram with all its internal lines on-shell with a MLT(1) diagram and one off-shell momentum. In the case of the N 2 MLT diagram, it is also decomposed into two terms. One of these terms is a convolution of an NMLT(3) diagram with a MLT(L − 3) diagram, and the other term is a convolution of one fully on-shell MLT(L − 3) diagram with a MLT(3) diagram with two external momenta inserted in the first and third internal set. These external momenta play a fundamental role in splitting the internal sets into two propagators each. In this manner, the following relations are fulfilled which justifies the factorisation relations presented in Sec. 2. The topological complexity of the subdiagrams in the convolutions is an additive quantity whose sum coincides with the topological complexity of the original diagram. Recalling that the family of Feynman diagrams with topological complexity k with L loops is N k−1 MLT(L), for the NMLT(L) has k NMLT = 2 and the term in the decomposition with no fully on-shell diagrams is given as the convolution of MLT diagrams (with k MLT = 1), thus 2k MLT = k NMLT . For the N 2 MLT(L) it is given k N 2 MLT = 3, while the term in its decomposition with no fully on-shell diagrams is the convolution of an NMLT(L) and a MLT(L) diagram, hence, These results suggest a factorisation formula for each topological complexity family. Indeed, these formulae can be reached through the nested residues and reveal a relation between one topological family with topologies with lower complexity. Within such an approach, every topological family can be expanded by factorising the MLT diagrams; thus, for any topological family, its causal structure can be expressed following Eq. (3.5). We profited from these structures to derive all-order formulae for N 3 MLT and N 4 MLT families, as explained in Sec. 6.

Asymptotic expansions within LTD
The LTD formalism is straightforwardly applicable to obtain several simplifications of the integrand-level representations of loop amplitudes, such as asymptotic expansions. The interest in asymptotic expansions arises from their potential to facilitate analytic results in specific kinematic configurations, particularly when full analytic calculations are not (yet) possible. An asymptotic expanded result can be of great interest since it showcases the relevant behaviour of the full amplitude in the needed kinematic limit. There are many observables where an analytic result is not necessary for every set of kinematics and where specific limits are the window to test potential discrepancies between experimental measurements and SM predictions, thus identifying new physics contributions. Furthermore, in the context of the local cancellation of IR singularities, expanded integrands could be very convenient to reduce computation time. While maintaining the correct analytic structure in the divergent limit and thus allowing for the combination with the real-emission contributions, the less complicated form of the expanded virtual contributions is expected to evaluate faster during the point-by-point process of numerical integration.
Since after the application of LTD an amplitude is reduced to an integral over Euclidean three-momenta, the size of the appearing scalar products can be unambiguously be compared to external scales. This provides a good starting point for the development of a well-defined formalism for asymptotic expansions of the integrand. Specific asymptotic expansions in the context of LTD have been presented for the first time for the process H → γγ at one loop [75]. Recently advances in the generalization of the formalism have been published in Refs. [76,77].
Due to the fact that scattering amplitudes are determined by their analytic properties general considerations for integrand-level expansions should start with the propagators as the origin of divergences. The numerator, while playing a role in the appearance of UVdivergences, is not relevant for the discussion of asymptotic expansions since within LTD the singular UV behaviour is neutralised through local renormalisation before integration.
We can reparametrise the dual propagators in the following form that is more suitable for asymptotic expansions If Γ ij +∆ ij vanishes, the dual propagator is not expanded. Otherwise the starting point for the asymptotic expansion is to demand that the condition be fulfilled for the whole range of the loop integration space except for potentially small regions around physical divergences. Note that since dual propagators only appear in integrands where one loop momentum has been set on shell, the condition has to be fulfilled in the Euclidean space of the loop three-momentum. The dual propagator can then be expanded as Further simplifications arise whenever k ji = 0. In that case, with the change of variables , the denominator of the expanded dual propagator takes the easily integrable form The intention to rewrite the denominator in this shape determines the parameters Γ ij and r ij appearing in the expansion to be restricted by the conditions assuming |r ij | ≤ 1. In the types of limits where one hard scale Q is available, Q 2 i = ±Q 2 can be identified. The sign is determined by the sign of the hard scale in the expression k 2 ji +m 2 i −m 2 j . This type of expansion facilitates the analytical integration based on integrals of the form In addition to the relations in Eq. (4.5) further conditions have to be respected by the expansion parameters to ensure that the expansion converges both at integrand-and at integral-level. Specifically, it is fundamental that the analytic behaviour of the dual propagator is not fundamentally changed, i.e. that for a propagator with a singularity the expansion also has to display that singularity, while the expansion of a non-singular propagator is to be finite throughout the whole integration domain. The infinitesimal imaginary prescription of r ij given in Eq. (4.5) accounts properly for the complex prescription of the original dual propagator and therefore its causal thresholds. A different approach has to be taken for threshold limits where a hard scale is not identifiable. In this case the pole position of the non-expanded propagator can be expanded to identify the correct r ij for the asymptotic expansion.
The formalism of expanding the dual propagator has been developed through its application on the locally renormalised scalar two-point function Applying Eq. (4.3) on the appearing dual propagators G D (q 1 ; ) and G D ( ; q 1 ) leads to a simplified expression that can be integrated without needing to specify what type of limit is considered, giving the very general result The coefficients c are simple functions of the appearing scales and are given in Ref. [77] just like the parameters needed for different limits. The expansion converges well both at integrand-and at integral-level in the limits of one large mass, a large external momentum as well as when approaching the physical threshold both from below and from above. Comparison with the established method of expansion by regions [103][104][105][106][107] has shown faster convergence as well as emphasised the advantage of decreasing degree in UV divergence with each order in the expansion. Renormalisation within our method thus only involves the lowest orders of the integrand-level expansion. Higher terms are optional for increasing precision and can be added straight-forwardly without the need to ensure further UV cancellation.
Using the same approach we have found integrand-level expansions for the scalar threepoint function (4.10) both for the limit of a large mass as well as for a small mass The Euclidean structure of the dual integrand, a( ), can be exploited into an even more direct way by applying Taylor expansions. Here it is important to use different assumptions on the size of the loop three-momentum in separate parts of the integration range. For the scalar two-point function in the large-mass limit this amounts to calculating The expanded integrand converges well in the full range of the loop momentum as can be seen in Fig. 7. Details on this method and on how to fix the matching scale lambda that separates the appearing so-called dual regions can be found in [77]. With the calculations above being performed in the traditional LTD representation we will finally show an example of using the manifestly causal representation to expand a multiloop integral. The MLT structure offers an ideal starting point, being void of unphysical singularities. For the large mass limit we find where s,0 . Notice that due to the lack of dependence on p 0 both in x L+1 nor in λ 0 L+1 the asymptotic integrals on the right-hand side of Eq. (4.14) are a function only of the internal masses, m s , to all loop orders. In this way, we are optimistic about future applications of the manifestly causal LTD representation from Refs. [92,[95][96][97] to speed up the calculation of efficient and smooth asymptotic expansions.

Manifestly causal representation and numerical efficiency
As explained in Secs. 2 and 3, the application of the nested residue strategy leads to manifestly causal integrand-level representations of any multi-loop multi-leg Feynman amplitude.
The main advantage of such a representation is the absence of non-physical singularities, because only terms compatible with causality remain in the final result.
Starting from the compact LTD representations of the NMLT and N 2 MLT multi-loop topologies presented in Ref. [92] in terms of nested residues, we reconstruct their full analytic expression in term of causal propagators only. This operation is performed by using numerical evaluation over finite fields [47,108], in which we use the C++ implementation of the FiniteFlow [109] algorithm together with its Mathematica interface. In the following Section, we show the explicit causal formulae for NMLT and N 2 MLT families for an arbitrary number of loops, and discuss their interpretation in terms of entangled causal thresholds. Numerical examples for 4-loop diagrams are presented, in order to provide a comparison with available results and test the efficiency of our approach.

Next-to-Maximal Loop Topology (NMLT)
At three-loops, the MLT family is not enough to describe the whole set of possible topologies. Thus we need to consider the NMLT and N 2 MLT topologies, whose general dual representations were explained in Sec. 2. To simplify the presentation, we start by considering NMLT configurations with one single propagator in each set and no external momenta. We need to add an additional internal line, whose momentum is given by A pictorial representation of NMLT is provided in Fig 1b. After computing the nested residues and adding all the contributions together, we get where the causal propagators are given by This expression was reconstructed by using numerical evaluations over finite fields, although partial fractioning leads to the same result within a similar computing time. In the following, we shall note that simplifications are not straightforward when dealing with more complicated topologies. An interesting consequence of the compact form of Eq. (5.2) is that it allows a reinterpretation in in terms of entangled causal thresholds. Each λ i is associated to a causal threshold singularity, which might take place when the momenta flows are oriented in the same direction. Thus, the product of causal denominators represents a configuration in which two (or more) sets of propagators can simultaneously go on shell. The factorisation of NMLT (and more complicated topologies) as products of MLT configurations is the reason behind this behaviour [92]. A graphical interpretation of the entangled causal structure of NMLT vacuum diagrams is provided in Fig. 8, with the dashed lines indicating the propagators that can be set simultaneously on shell. 5.2 Next-to-Next-to-Maximal Topology (N 2 MLT) As mentioned before, reaching a full description of N k MLT with k ≤ 2 is enough to obtain the causal representation of up to three-loop scattering amplitudes. The so-called N 2 MLT can be built recursively from the NMLT by adding an additional line with momentum The minimal example of such topology is the Mercedes-Benz diagram (L = 3) shown in Fig. 1c. For the sake of simplicity, we restrict here to the case without external momenta.
Using the LTD representation in Ref. [92], we can add together all the contributions and obtain  6) and N ({q i,0 . Compared to the NMLT case, the complexity of the polynomial in the denominator makes it highly non-trivial to unveil a formula similar to Eq. (5.2). Thus, we rely on the analytic reconstruction over finite fields to obtain It is worth appreciating that the package Lotty was used to efficiently reach these results [101]. As in the case of NMLT topologies, it is possible to interpret this result by using entangled thresholds. This time, there are products of three causal propagators. Notice that not all the combinations of causal propagators are allowed. This is because causal propagators exhibit some compatibility issues, which can be explained by digging into graph theory. More details about this issue and the connection with Cutkosky's rules can be found in Ref. [102].

Adding external momenta and higher-powers
For realistic multi-loop scattering amplitudes, we need to take into account external legs. The analytic reconstruction algorithm used to obtain compact formulae for vacuum diagrams can be also applied to topologies with external momenta. We want to highlight that the insertion of external momenta does not affect the causal physical behaviour of these integrals: the difference w.r.t. the vacuum case is that the entangled configurations are duplicated, according to the direction of the energy flow for external particles. For instance, a generic N 2 MLT with external legs inserted in the vertices is given by More details about the algorithms used to perform this computations can be found in Refs. [96,101].
Besides that, it was shown in Refs. [82,84,85,110] that the presence of self-energy insertions or generic scattering amplitudes, as well as some local UV counter-terms, might require to consider propagators with higher-powers. As explained in Refs. [96,97], the causal structure of these amplitudes can be obtained by applying a differential operator. Explicitly, we can raise the power of the propagators by taking derivatives w.r.t. q which suggests the definition of the operator The iterated application of this operator to the causal representation of scattering amplitudes with single powers of the denominators leads to causal representation of the corresponding multi-power amplitude, as carefully explained in Ref. [97].

Numerical implementations
Finally, we would like to highlight that the causal representation of multi-loop multi-leg Feynman integrals and scattering amplitudes leads to a very smooth numerical integration. In Ref. [96], we studied several examples and compared the performance of the LTDinspired approach with other standard frameworks (such as Fiesta 4.2 and SecDec 3.0). We considered the generic scalar amplitude

Universal opening at four-loops
The multi-loop topologies that appear for the first time at four loops are characterised by multi-loop diagrams with L + 4 and L + 5 sets of propagators which correspond to the next-to-next-to-next-to maximal loop topology (N 3 MLT) and next-to-next-to-next-tonext-to maximal loop topology (N 4 MLT). In fact, N 4 MLT naturally includes all N k−1 MLT configurations, with k ≤ 4.
The N 4 MLT family is represented with three main topologies, which were checked with QGRAF [111] and are shown in Fig. 11. Based in the similarity of these topologies with the insertion of a four-point subamplitude with trivalent vertices in a larger topology, a unified description is proposed. The three N 4 MLT topologies are interpreted as the t-, sand u-kinematic channels, respectively, of a universal topology.
The three topologies contain L + 4 common sets of propagators, and one extra set which is different for each of them. Each of the first L sets depends on one characteristic loop momentum s , the remaining four common sets are established as linear combinations of the loop momenta. The extra sets are the distinctive key to each of the channels in the universal topology where the momenta of their propagators is identified as different linear combinations of 2 , 3 and 4 .  To assemble the three N 4 MLT channels into a single topology, a current J is defined as the union of three sets that characterise each channel, J ≡ 23 ∪ 34 ∪ 24 . The dual opening of this topology fulfills a factorisation identity similar to those presented in Ref. [92] for NMLT and N 2 MLT also valid regardless of the internal configuration. This factorisation identity is depicted in Fig. 12 and is called the universal identity in view of the fact that it is the only master expression required to open to nondisjoint trees any scattering amplitude of up to four loops. It also enables to infer the causal structure of the complete topology by exploring the causal behaviour of its subtopologies. The term A where the bold symbol s is used to indicate that these contributions are those containing on-shell propagators in the J-sets. The first term on the r.h.s. of Eq. (6.4) and Eq. (6.5) consists of a four-loop N 2 MLT and three-loop NMLT subtopology respectively. These terms are describing dual trees where all the propagators with momenta in J remain off shell. The second term on the r.h.s. of Eq. (6.4) and Eq. (6.5) collects contributions that characterise the s, t or u channel shown in Fig. 11. To obtain these contributions a propagator in either set 23, 34 or 24 are set on shell.
The explicit expressions for the s, t and u channel arising from Eq. (6.4) and Eq. (6.5) are presented in Ref. [95] where all the results are consistent with the absence of disjoint trees. Repeated propagators from self-energy insertions were treated as single propagators raised to specific powers and are not considered to generate disjoint trees when the repeated propagator is set on shell.
Notice that the number of trees in the LTD forest can also be computed through the combinatorial exercise of selecting, from the full list of sets, all possible subsets of L elements that cannot generate disjoint trees. Nevertheless, the momentum flows of the onshell propagators can only be determined through the nested residues. Also, it is important to mention that the number of terms for a given N k−1 MLT topology in Eq. (6.3) scales with the number of loops and linearly with the number of propagators per loop set, but the sum over residues, equivalently over internal propagators, is implicitly accounted in this expression.

Causal representations
To confirm the causal conjecture for the N 4 MLT family, the strategy proposed in Ref. [96] is applied to the multi-loop N 3 MLT, t, s and u channels. Each topology considers a configuration with one internal propagator in each loop set, four external momenta for N 3 MLT and six external particles for t, s and u channels.
As a first step, the LTD representation is obtained for each of the selected topologies through the universal N 4 MLT expression in Eq. (6.3). After computing the nested residues and adding them all together a causal expression is found. The integrand of the causal dual representation reads in terms of on-shell energies, the energy components of the linear combination and causal denominators. Let us remember and emphasise that these causal denominators are constructed from sums of on-shell energies exclusively, and they represent potential singular configuration.
The results given by the straightforward application of LTD leads directly to a manifestly causal expression, however, the resulting numerator is a lengthy polynomial in the on-shell and external energies [95]. Therefore, a way to get a more suitable causal expression is to reinterpret it in terms of a number of entangled thresholds equal to the difference between the number of propagators and the number of loops by, e.g., analytical reconstruction from numerical evaluation over finite fields [47,108,109] as defined in Ref. [96].
The discussion was stated only for scalar integrals given that they fully encode all the compatible causal combinations. In general, tensor reduction commutes with LTD and can be used to deal with tensor integrals. Another feature exploited is that external momenta attached to interaction vertices that connect different loop sets do not alter the number of internal propagators and therefore the complexity of the causal representation. Therefore, the full causal expressions with external momenta can be deduced from the causal representation of the vacuum configuration by matching the momentum flows of the entangled thresholds.
Starting with the multi-loop N 3 MLT, there are thirteen causal denominators which are depicted in Fig. 13. The analytically reconstruction was done by matching all combinations of four thresholds that are causally compatible to each other.
Going forward to the causal representation of the N 4 MLT family, we have to consider all the entangled configurations with the presence of five causal thresholds. The t channel depends on the causal denominators already defined for the N 3 MLT configuration and additional nine extra causal denominators that depend on q (+) 23,0 where the corresponding configurations are shown in Fig. 14.
The s-channel can be obtained through the structure of the t-channel. The causal denominators are obtained just by a clockwise rotation of the t-channel and therefore by a permutation of the arguments of the causal denominators that are channel specific. Similar to the s-channel, the u-channel is obtained from the t-channel through the substitution 23 → 24 and by the exchange 3 ↔ 4 or 2 ↔ 234 (123 remains invariant). There are, however, three new configurations that arise because the u-channel is nonplanar. These new configurations are shown in Fig. 15.
One advantage of the causal representation is that the number of terms for a given N k−1 MLT topology is independent of the number of loops in the causal representation but requires to specify additional causal thresholds and additional causal entanglements when more internal propagators are considered. Furthermore, the most significant advantage of the causal representation with respect to the LTD representation stems from the core difference between them, the presence or absence of noncausal singularities. The straightforward application of the nested residue generates multiple threshold singularities, nevertheless, with a clever analytical rearrangement, the absence of noncausal singularities is achieved and leads to a causal representation which is more efficient and stable numerically in all the integration domain [95].

Novel developments on causality
The applications of the LTD-inspired methods have been spreading very fast in the last years. In particular, it turned out to be an excellent tool to unveil the structure of causal singularities of multi-loop multi-leg Feynman integrals and scattering amplitudes. There have been very recent findings regarding general all-loop formulae to describe any N k MLT amplitude, by using clever algebraic relations among them [100]. This computational technology has been implemented in the package Lotty [101], which allows to automatically obtain the causal representation of multi-loop multi-loop Feynman integrals and scattering amplitudes. On the other hand, investigations inquiring on the geometrical aspects of multi-loop diagrams were performed. In Ref. [102], we used concepts from graph theory to describe all the possible causal propagators and the allowed entangled thresholds associated to a given diagram. Describing diagrams in terms of cusps (i.e. interaction vertices) and edges (i.e. sets of propagators connecting two fixed vertices) was also proposed in Ref. [100]. This description leads to the cusp matrix, a geometrical object that allows one to completely characterise the causal structure of a given diagram. Explicitly, we realised that the causal propagators are associated to all the possible connected binary partitions of cusps. This has a very nice connection to a generalised geometrical interpretation of the Cutkosky's rules [88]. Also, we implemented a computational algorithm that benefits from graph theory tools to efficiently obtain all the causal propagators. We tested it against the direct nested residue calculation and subsequent identification of λ's: the geometric algorithm was, at least, a factor ten faster.
The information contained in the cusp matrix allow us to exactly reconstruct all the possible entangled thresholds by imposing geometrical selection criteria. These can be summarised in the following list: 1. All the lines are cut: Each possible entangled combination of causal denominators must involve the on-shell energies of all the propagators. i,0 ) and the energy component of the external momenta is determined by the orientation matrix, which can be built from the cusp matrix.
In Fig. 16, we show examples of causal propagators which can not be entangled for the Mercedes-Benz diagram (i.e. a particular case of N 2 MLT or a four-cusp topology). On the left, the criteria 2 is not fulfilled, since the associated partitions of cusps are not disjoint. On the right, we consider three causal propagators which can not be consistently oriented (criteria 3): the edge connecting the cusps 1 and 3 does not allow a compatible orientation.
Thus, given a diagram, determining its cusp matrix and imposing the selection criteria 1-4, we can obtain a formula describing its causal structure. In accordance with the results presented in Ref. [100], we found that describes the causal structure of any multi-loop multi-leg Feynman diagram or collection of Feynman diagrams with topological complexity k − 1. The set Σ contains all the subsets of products of k causal propagators λ ± i fulfilling the selection criteria 1-4, and N σ is given by the application of an operator depending on the subset of σ. Thus, this result confirms that multi-loop multi-leg amplitudes can be purely described in terms of cusps and edges, regardless the number of loops or external lines.

Four-Dimensional Unsubtration (FDU)
The standard computation of accurate predictions at Next-to-Leading Order (NLO) in QFT requires to deal with non trivial integrands. These integrands carry multiple singularities due to the low (IR) and high (UV) energy regimes. The cancellation of those divergent quantities, typically occurs at integral-level. By considering the Kinoshita-Lee-Nauenberg (KLN) theorem [112,113], soft and collinear divergences (IR singularities) are removed, and by adding suitable ultraviolet counter-terms, UV singularities are renormalised, rendering the cross section finite. This algorithm for cancelling singularities can be extended to Next-to-Next-to-Leading order (NNLO) and beyond; however, since the complexity of the integrands is higher at higher orders, this procedure is reaching a bottleneck because cancellation of singularities must be achieved at integral level, i.e., the → 0 limit is taken once the cross-section is known in DREG. Since LTD transform loop integrals in phase-space integrals, among other properties discussed in this document, the cancellation of singularities could occurs at integrand level. Hence, based on the LTD theorem, the Four Dimensional Unsubtraction [82][83][84][85]110] (FDU) method presents a new paradigm to compute observables in four dimensions, since the definition of cross-section is free of IR and UV singularities. The first application of the LTD at cross-section level was done in a toy model based on the simplest φ 3 theory. In Ref. [82] it was shown for the first time that LTD is powerful to build physical observables in four dimensions by a proper mapping between real and virtual kinematical variables.

Local cancellation of infrared singularities within FDU
Let us start the discussion of the FDU approach by analysing the cancellation of soft and collinear divergences. The simplest scenario is the decay process 1 → n where the Born level cross-section is given by, with |M (0) n | 2 the LO contribution and S 0 is the IR-safe measure function. The virtual correction to the LO cross-section is computed as where M (0) n |M (1) n represents the interference between the Born level and one-loop amplitudes. It is important to emphasise at this point that FDU requires also to maintain selfenergy contributions because they contain both IR and UV divergences that will contribute to the full cancellation of singularities. The IR singularities from the virtual component are cancelled against real radiation contributions.
For the sake of simplicity, let us consider that the process under study only has finalstate radiation singularities. Hence, the real contribution given by contains all Feynman diagrams with one extra particle in the final state in M (0) n+1 and S 1 is the measure function for n + 1 particles. It is important to recall at this moment that primed variables are used to describe the momenta of real matrix elements while unprimed variables are labelling Born and loop level momenta.
The FDU starts with the analysis of the virtual cross-section which can be decomposed as a sum of dual contributions as, with I i (q i ) the N dual integrals that arise from the direct application of the LTD theorem. We notice that each I i (q i ) has set one internal particle on shell and it is characterised by q i . Then, the loop integral has been transformed to phase-space integrals, therefore each dual integral behaves as the contribution of one real particle emitted. In order to achieve the complete cancellation of IR singularities, each σ D,i must be paired, at integrand level with a corresponding real radiation twin component. The real radiation is sliced as, where R i represents a partition of the full integral, such that i σ (1) R and considering that only one IR divergent configuration is allowed in each R i . Finally, in order to merge both integrals, a mapping between the kinematical variables is implemented. On the one hand, each partition in σ (1) D,i contains n external momenta plus one extra integration variable corresponding to q i , therefore, each σ (1) D,i is mapped into its real sliced contribution as, where T i is a bijective transformation among real and virtual variables. The momentum mapping between variables is analogous to the one used by the dipole method [57,58], since the singularities are associated to soft emissions or the double collinear limit [114,115]. These singular IR configurations are associated to specific contributions in the virtual part, by selecting two massless partons per partition. On one hand the spectator and on the other hand the emitter. Hence, the four momenta of the emitter and spectator in addition with the loop three-momentum are used to reconstruct the kinematical phase-space of the real emission cross-section where there is a similar configuration, i.e., an emitter decaying into two partons in a soft or collinear regime. Explicitly, if p i is the momentum of the final-state emitter, q i the internal on-shell momentum prior the emitter and p j is the momentum of the final-state spectator, we apply the momentum mapping with p r the momentum of the extra radiation of the process and, in the massless scenario, . Furthermore, this mapping preserves momentum conservation since p i + p j + k =i,j p k = p i + p j + p r + k =i,j p k is fulfilled. A similar mapping is used when considering massive particles and it was shown in Ref [85]. The extension of this formalism to higher-orders (i.e. NNLO and beyond) will require additional kinematical transformations which take care of the singular behaviour of scattering amplitudes in the multiple-collinear limit [116][117][118][119][120][121][122].

Self-energy insertions and renormalisation
Before moving forward with IR cancellation, let us review the local renormalisation of UV singularities. The standard cancellation of UV singularities requires the renormalisation of field wave-functions and couplings. In the FDU formalism, this feature is obtained through the construction of local UV counter-terms. At one-loop, the scenario where massless and massive particles are propagating in the loop has been studied [84,85]. Remarkably, it was shown that a smooth transition between massless and massive renormalisation constants takes place and UV singularities are well understood. Let us highlight a crucial difference between the standard renormalisation constant in DREG and LTD. Wave-function renormalisation constants are obtained from self-energy diagrams. In particular, massless bubble diagrams in DREG are neglected in the renormalisation procedure since IR and UV divergences are considered as equal and they cancel out. However, the same divergent poles in the FDU formalism contribute separately. Specifically, it means that IR singularities of loop diagrams vanish with only the IR poles of the real radiation and there are no mixed cancellation between UV and IR poles. Therefore, the remaining UV divergences must be removed when renormalisation is implemented in the FDU scheme. Hence, we stress that integrands in the FDU are separated into the IR and UV domains, and this identification is crucial to render cross-sections free of singularities in the FDU formalism. We analyze the construction of the renormalisation constants in the FDU framework. By using standard Feynman rules, the unintegrated massive wave-function renormalisation constant, in the Feynman gauge, at one-loop is given by, Eq. (8.8) is the most general wave-function renormalisation constant since it includes the massless and massive case and, as it was previously mentioned, the transition to the massless case, ∆Z 2 (p, 0), is straightforward. The function ∆Z 2 (p 1 ; M ) contains singularities associated to the UV domain, therefore we must find the UV component of Eq. (8.8), ∆Z UV 2 , and subtract it in order to find a UV-free wave-function renormalisation constant, ∆Z IR 2 . The UV part is extracted by performing an expansion of the integrand around the UV propagator G F (q UV ) = (q 2 UV − µ 2 UV + ı0) −1 . In particular, for Eq. (8.8), it is found Hence, we define ∆Z IR 2 as which is free of UV singularities and the IR singularities poles are those needed to remove the remaining singularities from the virtual and real mapping of the cross-section.
To remove UV singularities it is important to build proper UV counter-terms that can be extracted by the direct study of the UV properties of physical amplitudes. After the combination of real and virtual matrix elements and renormalisation, the remaining amplitude has only UV singularities, in the following, we consider generic amplitudes with only UV singularities. At two loops, a generic two loop amplitude can be written as, where the integrand is a function of the loop variables 1 and 2 . UV divergences shall appear when the three-momenta | 1 | and | 2 | tend to infinity. In the two-loop scenario, there are three possible UV limits: i) | 1 | → ∞ and | 2 | remains fixed, ii) | 2 | → ∞ and | 1 | is fixed and, iii) | 1,2 | → ∞. In order to extract the UV behaviour of the integrand, we impose the replacement, for a given loop momentum j and then, we expand the expression up to logarithmic order around the UV propagator. This action is represented by the L λ operator. Notice, however, that the result shall generate a finite part after integration that has to be fixed to find the right value of the integral. Therefore, the first counter-terms is computed by with d j,UV the fixing parameter which makes the finite part of integral to be zero in the MS scheme. The remaining divergences shall occur when both | 1 | and | 2 | approach to infinity simultaneously. Then, to build a counter-term that mimics this behaviour, the following replacement is implemented, Figure 17. Momentum configuration of NLO QCD corrections to the process A * → qq(q), assuming that the decaying particle does not couple to gluons.
on the subtracted integrand. Similarly to the previous UV counter-term, the application of the L λ operation shall produce a finite piece that has to be fixed to build properly the counter-term, A with d UV 2 the fixing parameter of the double UV limit. Having all UV divergences under control, the renormalised amplitude at two-loops, A R , can be constructed by the subtraction of all UV counter-terms, such that is free of IR and UV singularities.
In the multi-loop scenario, multiple ultraviolet poles will appear, since all loops could tend to infinity at different speed.

A * → qq at NLO in QCD
The FDU method has been studied in different processes, for the sake of simplicity, we recall the application of the algorithm in A * → qq a NLO in QCD, with A = H, γ, Z.
At NLO, virtual correction A * (p) → q(p 1 ) +q(p 2 ) must be mapped into the real radiation process A * (p) → q(p 1 ) +q(p 2 ) + g(p r ). As it was explained in previous subsections, we must also keep all virtual contributions since self-energy contributions are crucial to cancel out all IR singularities. In Fig. 17, we present the complete set of Feynman diagrams of the real and virtual processes.
It is important to recall that in this example we consider massive quarks, i.e., p 2 1 = p 2 2 = M 2 and p 2 1 = p 2 2 = M 2 and p 2 r = 0. It is convenient to write the momenta p 1 and p 2 in terms of massless four vectors,p 1 andp 2 , such as with β ± = (1 ± β)/2, β = √ 1 − m 2 and m = 2M/ √ s 12 . Therefore, since the momentum of the gluon is labelled in the virtual contribution as q 1 and LTD shall set it on-shell, q 2 1 = 0, the momentum mapping proposed for the merging of virtual and real cross-sections is where α 1 and γ 1 are parameters determined by imposing on-shellness. Notice that Eq. (8.18) fulfills momentum conservation. This mapping will be used to remove some IR singularities of the integral related toδ(q 1 ), the remaining IR singularities shall cancel through renormalisation. Now, since there are other singular IR behaviours in the real that can be recognised in the virtual, when q 2 2 = M 2 , the second momentum mapping is proposed to be, with α 2 and γ 2 parameters which are also determined by the on-shell conditions. In virtue of soft and collinear divergences in the real cross-section being only in this two LTD cuts, we proceed to slice the full real element with the general function, 20) which in this particular case are explicitly, θ(y 2r < y 1r ) + θ(y 1r < y 2r ) = 1 , with y ir = 2p i ·p r /s 12 . The integration domain is well determined and it is shown in Fig. 18. We remark that this arrangement of the real element phase-space is general and it can be applied to any 1 → 3 real cross-section that is mapped into the 1 → 2 virtual cross-section.
Let us now compute the virtual decay rate of A * → qq. The LTD representation is given by, with the renormalised one loop amplitude computed as, and |M (1,UV) A is the unintegrated UV counter-term of the one-loop vertex correction. The real element, after the splitting of the phase-space in two domains, is written as, Figure 18. The dual integration regions in the loop three-momentum space, with ξ ⊥ = ξ 2 x + ξ 2 y . A , are presented in Fig. 19. The agreement with the analytical predictions is excellent in all cases. Furthermore, we notice that the massless quark scenario is recovered smoothly in the FDU however, in DREG this is not the case since individual contributions are not smoothly defined in that limit.

Local dual counter-terms for cross-sections
In the following, we will focus on the NLO QCD corrections to cross-sections involving massless colored particles and any number of other colourless SM particles. The extension of the following ideas to the massive case has been worked out following the same reasoning, which is presented in Ref. [86]. The reduction of any one-loop amplitude to a combination of scalar integrals is an easy task compared to the two-loop case. Of course, such a reduction does not alter the structure of the divergences of the amplitude but to properly reproduce all the IR singularities, in addition to the wave function renormalisation, one has to keep also the massless bubble integrals that are formally discarded usually (as already explained in Sec. 8). Once all the box functions are expressed in terms of their corresponding 6−dimensional version, it turns out that the IR divergences of a renormalised one-loop amplitude interfered with the leading order are given by the triangle scalar functions with two external partons connected by a massless propagator times the corresponding color connected tree level amplitude multiplied s ij = 2 p i ·p j , the massless bubbles times (1/2) the Born for every external (gluon) quark external particle multiplied times the corresponding Casimir, and the external wave function renormalisation constant time the Born. The above set of IR singularities is universal and, in fact, it reproduces the known formula for the poles of any one-loop amplitude in QCD. Inspired by the LTD theorem, we considered the cuts of these loop integrals that capture their infrared and collinear divergences. As a result of this procedure one obtains a novel subtraction scheme. For the sake of simplicity we will restrict the discussion to the case of a colourless initial state and so we will consider colored radiation exchanged only among the final state partons. Our general formula for the dual counter-term in the case of a quark emitting a gluon is the following: 1) where N in represents all the non-QCD factors, including the 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 ac and b, respectively. The functions V qg,b and G qg,b in Eq.(9.1) are defined by where V qg,b groups together the dual contributions from the reduction of the virtual amplitude and G qg,b is the contribution of the quark wave function renormalisation. In the above formulas p i and p j represent the momentum of the emitter (with color index ac) and spectator (color index b) in the Born kinematic respectively, while q 1 is the loop momentum that is associated to the gluon in the real radiation kinematic. In the case a gluon is radiated collinear to another gluon the corresponding dual counter-term σ DS gg,b is given by where the symmetrisation p a ↔ p c has to be performed on the right-hand side. Once again, in Eq. (9.3) we have extracted only the dual contributions with q 1 on-shell, obtaining Finally, the case of a gluon splitting into a collinear quark-antiquark pair receive contributions only from the wave function contribution so that the corresponding dual counter-term σ DS qq,b is given by (9.6) Once a proper momentum mapping is considered, these counter-terms can be analytically integrated over the region in which p a · p c < p b · p c in the real phase space and added back to the virtual contribution. This last operation corresponds to extract all the infrared and collinear divergences out of the virtual matrix element. We have used the mapping adopted in the seminal work by Catani and Seymour [58].

γ * → 3 jets at NLO
As an example of application, we consider the NLO correction to the three-jet production in e + e − annihilation γ * → q(p 1 )q(p 2 ) g(p 3 ). For simplicity we will restrict to the case of gluon radiation, namely γ * → q(p 1 )q(p 2 ) g(p 3 ) g(p 4 ). The case of four quark production can be treated along the same lines. 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. Once a real kinematic configuration (p 1 , p 2 , p 3 , p 4 ) is generated, the six kinematic invariants s ij are analyzed and as a consequence six dual counter-terms are activated. To each emitter-spectator pair, we associate the following dual counter-terms:  Since the process involves only three colored partons, the colour algebra factorises and one has T 1 T 2 = C A /2 − C F and T 1 T 3 = T 2 T 3 = −C A /2. In the virtual sector, we use the integrated version of the dual counter-terms that are given by e + e− annihilation at √ s = 125 GeV. Jets are clustered using the Durham jet algorithm with resolution parameter set at y cut = 0.05. A perfect agreement with the same computation performed using Catani-Seymour dipoles is observed. The same level of agreement is observed for any other differential distribution.
In this way, we have shown that a fully local implementation of an LTD-inspired method for computing cross-sections is reliable. These studies complement the ones presented in Sec. 8, indicating that the FDU framework constitutes a powerful strategy for reaching higher-order corrections with fully numerical methods. An extension of this technology to calculate N 2 LO corrections for IR-safe observables is being actively studied.

Outlook and further developments
In this manuscript, we reported on the recent developments in the numerical evaluation of multi-loop scattering amplitudes through the application of the Loop-Tree Duality (LTD) formalism. Based on the original approach, at one and two loops [68,70], we presented a novel formulation of LTD [92], which allowed to bootstrap the treatment of Feynman integrals, regardless of the loop order, and shed light to a complete automation. In effect, the decomposition of two-and three-loop scattering amplitudes in terms of dual integrands was obtained by introducing the concepts of maximal (MLT), next-to-maximal (NMLT) and next-to-next-to-maximal (N 2 MLT) loop topologies. These families of diagrams share general properties that allow allows for the generalization of the dual representation at any loop order.
Since LTD heavily relies on the application of the Cauchy residue theorem, we studied in details the treatment of multivariate rational functions in Ref. [97], giving, in this way, a mathematical support to the novel formulation of Ref. [92]. Likewise, we proved all conjectures that were provided in the latter, thus, rendering in a well organised formalism ready to focus on physical applications. Hence, in view of the LTD decompositions of up-to three-loop scattering amplitudes, the natural extension, along the lines of Ref. [92], was considering the four-loop case [95]. In effect, the complexity in the treatment of four-loop topologies increased because, differently from the two-and three-loop cases, this loop order was not described by only one loop topology, as in the lower cases. Hence, we introduced new configurations that appears at four loops: the N 3 MLT and N 4 MLT families. In this way, our formalism allowed a complete understanding and decomposition of the dual representation of any up-to four-loop scattering amplitude.
In the spirit of profiting from the dual representation of integrands through LTD, we observed that the complete sum of the latter, for a given Feynman integral, led to a representation of an integrand that only manifests physical information. In other words, unphysical singularities or pseudo-thresholds cancel out when all dual integrands are summed up. This is a remarkable feature of LTD that was originally observed in Refs. [79,87,92] and extensively studied in Ref. [96] for MLT, NMLT and N 2 MLT configurations. In fact, the work performed in Ref. [96] made a comparison between the numerical stability of dual and causal integrands. While in the former numerical instabilities, originated from pseudo-thresholds, were present, the latter displayed a smooth behaviour: the numerical evaluation was much more stable and this allowed a straightforward computation of the Feynman integrals. In particular, we considered the numerical integration of ultraviolet finite integrals at three and four loops with presence of several kinematic variables.
The causal representation of scattering amplitudes has recently been a very active topic due to the simple structure of these integrands in terms of causal (physical) thresholds. In fact, an understanding of the structure of the latter started to appear in the literature by means of geometric properties [102] and is yet under consideration. On top of the former approach, it was recently conjectured a breathtaking formulation of the all-loop causal representation of multi-loop scalar integrands, obtained from the features that characterised a topology, cusps and edges [100]. This formulation provides, regardless of the loop order, the most symmetric causal representation of Feynman integrands, in a complementary and independent way to the nested residue calculation used by the LTD formalism. Thus, to give a support to these conjectures, we provide the Mathematica package Lotty to automate both dual and causal representation of scattering amplitudes [101].
From the above-mentioned discussion, we claim that a treatment of Feynman integrands at an arbitrary number of loops is under control. Moreover, this is not the end of the story when comparing theoretical predictions with data derived by collider experiments. In effect, one yet has to deal with ultraviolet (UV) and infrared (IR) singularities. While the framework at next-to-leading order (NLO) is completely understood by means of the Four-Dimensional Unsubtraction scheme [82,84,85] and the dual subtractions [86], a clear path towards calculations at higher orders (NNLO and beyond) needs to be devised. Preliminary studies of IR-safe scattering amplitudes has been recently considered in Refs. [79,81]. Besides, the extension to the calculation of cross sections at NNLO is under study, in particular, focusing on a careful local treatment of IR singularities.
We expect that the ideas presented in this manuscript will allow to elucidate in detail the state-of-the-art of the LTD formalism and the causal representation of scattering amplitudes by means of the latter. Likewise, the ideas that we have developed will certainly supply the scientific community with a powerful strategy to tackle the calculation of multiloop Feynman integrals, in which phenomenological and formal applications can efficiently be carried out.