Loop Tree Duality for multi-loop numerical integration

Loop Tree Duality (LTD) offers a promising avenue to numerically integrate multi-loop integrals directly in momentum space. It is well-established at one loop, but there have been only sparse numerical results at two loops. We provide a formal derivation for a novel multi-loop LTD expression and study its threshold singularity structure. We apply our findings numerically to a diverse set of up to four-loop finite topologies with kinematics for which no contour deformation is needed. We also lay down the ground work for constructing such a deformation. Our results serve as an important stepping stone towards a generalised and efficient numerical implementation of LTD, applicable to the computation of virtual corrections.


I. INTRODUCTION
Loop integrals are an essential component of fixedorder corrections to collider cross-sections. Analytic techniques enjoyed a durable success in this matter, but it has becomes increasingly evident that a further breakthrough necessitates a radical change of perspective. Numerical approaches are a promising alternative and have already been extensively explored for Feynman amplitudes using sector decomposition (see e.g. [1][2][3][4][5][6][7]). More recently, direct integration of finite loop integrals in fourdimensional Minkowskian momentum space have been considered, together with the necessary complex contour deformation for handling integrable threshold singularities [8][9][10][11]. In this letter we study the possibility of rewriting an n−loop integral as a sum of terms with n additional on-shell conditions by analytically integrating over loop energies using residue theorem. The ensuing identity is called Loop Tree Duality [12] (LTD). LTD is appealing from a numerical standpoint for at least four reasons: (1) the n−loop integral dimensionality is fixed to 3n irrespective of the topology considered, (2) integrable singularities can be shown to be confined to a bounded volume [13] and are absent when considering certain kinematical configurations, (3) momentum-space divergent integrals naturally lend themselves to be regularised with local UV and IR counterterms [14][15][16][17][18][19][20][21][22][23] or even (4) through a direct combination with the corresponding real-emission contributions in the case of physical amplitudes [24][25][26].
In this work we derive a novel multi-loop LTD expression by iteratively applying the residue theorem, carefully keeping track of the propagation of Feynman's causal prescription. This differs from the expression derived in ref. [27] where distributional identities between Feynman and dual propagators are used. It also differs from the work of ref. [28] which seeks to achieve a similar goal using a multidimensional version of the residue theorem, arriving at an expression which is incompatible with ours.
We proceed to numerically apply our LTD construction to various scalar loop topologies ranging from one to four loops. In all cases, we find agreement to better than 1%, thereby validating our procedure and the expected structure of integrable as well as cancelling singularities exhibited by each term of our LTD expression. We also determine the constraints on a contour deformation and construct deformation vectors satisfying them. A first preliminary result is given for a two-loop LTD integration using a contour deformation.
The outline of this work is as follows. In section II we derive our LTD expression for an arbitrary number of loops. In section III we discuss our numerical implementation and results. In section IV we present first steps towards a general contour deformation. Finally, we give our conclusions in section V.

II. LOOP TREE DUALITY FORMALISM
We consider the following general expression for an nloop integral where e is the set of indices labelling the edges of a Feynman diagram and the numerator N is a regular function of the loop momenta. The Feynman propagator 1/D i depends on the four-momentum q i ≡ (q 0 i , q i ), the mass m i and the positive causal prescription iδ. We consider pairwise distinct Feynman propagators, each with two first order poles in q 0 i located at σE i ≡ σ q 2 i + m 2 i − iδ, σ ∈ {±1}, where +E i lies in lower complex half-plane. Introducing the signature vector s i = (s i1 , . . . , s in ), s ij ∈ {±1, 0}, we write q µ i = n j=1 s ij k µ j + p µ i , where p µ i is a shift that depends on external momenta.
The integration over the momenta k j can be split up in an integration over a spatial part k j and the energy k 0 j . We now derive our LTD formula by per-forming the energy integrations one after the other, following an arbitrary fixed order of the energy variables k 0 = (k 0 1 , . . . , k 0 n ). We construct this iterative procedure by considering a contour for each energy integration variable along the real line, and closing on an arc in either the upper (with winding number Γ j = +1) or the lower (Γ j = −1) complex half-plane. We assume the integral along the arc to vanish, such that the integral along the real line equals the sum of residues at poles located within the contour. The iterative computation of each loop energy integration yields where we introduce the set of ordered lists of edge indices I = {(i 1 , . . . , i n ) ∈ e n | det (s ij j ) 1≤j≤r = 0, ∀r ≤ n}, which guarantees that for every iteration j, where we integrate out k 0 j , the propagator labelled by i j depends on k 0 j . Note that this set can contain several permutations of the same indices. The residue of f = N/ i∈e D i is being evaluated at the pole locations implicitly defined through the solutions k 0 = k σ i to the following linear system where the signature matrix s i is a totally unimodular matrix. Each residue contributes to the integral if the pole location is within all contours of energy integrals already performed, corresponding to the condition Γ r [k σ i,r ] > 0, ∀r ≤ n. The imaginary part of the poles in the energy variable k 0 r is computed using Cramer's rule for the last row of the subsystem of (4) arising after every iteration. Its final expression is given by which explicitly shows how the imaginary parts of poles selected by previous iterations propagate to the imaginary part of the pole contributing at iteration r. Eq. (2) contains Heaviside functions Θ with complicated arguments. However, we have checked that for all topologies from one to six loops the Heaviside functions that do not identically evaluate to either 0 or 1 cancel pairwise. In fact, we find that for each loop momentum basis of the corresponding loop graph, only one combination of energy signs contributes to I with a definite prefactor (−i) n . We call this combination of signs the cut structure. Therefore, we conjecture that eq. (1) can be written as where b is the set of all edge indices labelling a loop momentum basis and B is the set of these sets for all loop momentum bases. The cut structure of b is denoted with σ b . Each loop momentum basis is assigned a residue, henceforth referred to as a dual integrand, reading where solving {q 0 Note that the dual integrand is invariant under permutations of the elements in b, unlike eq. (3) that depends on the ordering within i. Furthermore, the complement t = e\b is the spanning tree of the graph. There is a one-to-one correspondence between a spanning tree t and a loop momentum basis b, hence the name Loop Tree Duality.
We expect the sum of residues obtained from analytic integration of loop energies to be independent of the specific loop momentum routing as well as choice of contour closure for each loop energy integration. We verified that these expectations are met by explicitly applying eq. 6 for various choices of routing and contour closures, each time retrieving the same numerical result for the sum of residues b∈B Res b [f ] for given numerical inputs k j . In performing these checks, it was convenient to have the cut structure construction algorithm automated and we provide the corresponding Python implementation as ancillary material.
and still features singularities if it can go on-shell. The inverse dual propagator vanishes on two singular surfaces where σ ∈ {±1}, and where s b ij and p 0,b i are defined implicitly through the change of basis q 0 The singular surfaces can be separated into two classes, which we call E-and H-surfaces. To distinguish them, we define the surface signs for the surface A singular surface where all surface signs are equal is called an Esurface, since its defining equation is the one of an ellipsoid when n − 1 loop momenta are kept fixed. Otherwise, it is called an H-surface, since its equation is the one of a hyperboloid when viewed as a function of at least one loop momentum.
We now provide the multi-loop existence conditions for H-surfaces. In the one-loop case, and in general when and the H-surface exists for real masses and loop momenta iff as already found in ref. [29]. In the case of |S b,σ i | > 2 and if exactly one H-surface sign differs from the others, whose index in b ∪ {i} we labelẽ, we define the following quantity: and the corresponding H-surface exists iff or when the surface signs contain at least two positive and at least two negative members.
The singularities of dual integrands on H-surfaces cancel pairwise in their sum b∈B Res b [f ], due to a mechanism referred to as dual cancellations [13,29,30], independently of the regulator δ. We checked both numerically and analytically that eq. 6 maintains the dual cancellation pattern of H-surfaces also beyond two loops.
E-surfaces satisfy the following existence conditions for real masses and loop momenta: We note that when the bound above is saturated, the E-surface is said to be pinched and it corresponds to the location of physical soft and collinear singularities of the loop integral which would require dedicated local counterterms for its regularisation. The singularities on existing E-surfaces must be regularised through a contour deformation satisfying its corresponding δ-prescription. We derive this prescription by writing the leading term of the Taylor expansion in δ of the imaginary part of ∆ σ,b i : We note that for E-surfaces we have the definite sign sgn [∆ σ,b i ] = −σ independent of loop kinematics. If no E-surface existence condition is satisfied, the integrand b∈B Res b [f ] has no singularities and it is therefore independent of the regulator δ. In this case, the numerical integration can be performed without a contour deformation, a feature that has already been shown at one loop in ref. [29] and two loops in ref. [31]. A first preliminary result for a two-loop LTD integration using a contour deformation will be given in section IV.

B. Discussion of previous work
In ref. [27] an alternative multi-loop LTD expression is derived by using distributional identities between dual and Feynman propagators instead of applying residue theorem. The main distinction of the final expression lies in the iδ-prescription of the dual propagator when diagrams beyond one loop are considered: in the case of ref. [27], the dual prescription is not equivalent to evaluating the on-shell conditions with complex energies q 2 i + m 2 − iδ and it is not a quantity independent of the order of integration in the loop variables, unless dual integrands with more on-shell conditions than loops are added.
The careful propagation of Feynman's causal prescription in the iterative approach discussed in the previous section is instrumental for obtaining a correct LTD expression for n−loop integrals. In the work of ref. [28], an alternative LTD construction is presented, where an averaging procedure over all contour closures is considered, invoking the multidimensional residue theorem. The imaginary parts of each propagator are taken to be independent of each other throughout the induction proof, thereby not considering the interplay stemming from taking multiple on-shell conditions, such as the one reflected in eq. 5.
Our construction allows for arbitrarily choosing to close the contour of each energy variable either in the upper or lower complex half-plane. Using this, we have explicitly constructed the expression resulting from averaging over all possible contour choices and we find a different combination of residues than the one reported in ref. [28]. It thus appears that the aforementioned interplay in the determination of the sign of the imaginary part of each pole does not disappear upon this averaging procedure. The integrand stemming from the direct combination of the integrands for each spanning tree given in ref. [28] does not reproduce our sum of dual integrands. Furthermore, we observe that this combined integrand does not realise dual cancellations. We therefore conclude that the LTD expression presented in ref. [28] is incorrect beyond one loop, to the best of our understanding.

III. NUMERICAL APPLICATION
LTD has shown to yield promising results at one loop [29] and has the advantage of not necessitating any computationally demanding symbolic treatment of the integrand and/or its numerator. This is different from sector decomposition techniques, which require building the Feynman representation of loop integrals together with the identification of sectors. Moreover, integration in momentum space is particularly appealing for its optimal scaling with the number of contributing scales. Compared to the 4d momentum space integration method described in ref. [9], LTD has at least 5 advantages: (1) the dimension of the integration is reduced to 3 per loop, (2) a complex contour deformation only needs to be applied on bounded E-surfaces, (3) masses do not complicate the contour deformation much, (4) specific kinematical configurations can be integrated without any deformation, and (5) its singularity structure can directly be related to real-emission contributions [32].
In this work we are mostly interested in demonstrating LTD viability for numerical multi-loop computations and in assessing the validity of eq. 6. Therefore, we apply LTD to loop integrals with external kinematics that do not yield singular E-surfaces, such that no complex contour deformation is required. This scenario offers a reliable numerical check of our LTD cut structures and of the numerical stability of the dual cancellations. Our implementation is a first important step towards handling loop integrals in the physical regime, which we briefly discuss in section IV.    We selected eight very different loop topologies, displayed in table II, to showcase the generality of the method. We report our results in table I and figure 1, with additional information (such as the exact input kinematics) given as ancillary material to ensure reproducibility of our work. The reference results are taken from the analytic expression for the four-point integrals [34], from Forcer [35] for two-point integrals, from Mad-Loop [33,36] for the decagon and triacontagon and Py-SecDec [7] for the six-and eight-point integrals (in which case the numerical error is also reported). We find perfect agreement in all cases, but note that scalar integrals whose superficial degree of UV divergence is -2 (II.f and II.h) are numerically more challenging. This is made manifest for example when comparing LTD results obtained for the loops II.f and II.g. We find no notable sensitivity of the numerical convergence to the external momenta multiplicity, internal masses or non-planarity of the loop graph.
For all eight benchmark loop integrals, we have explicitly verified that dual cancellations hold by sampling points on the H-surfaces for which we found that the sum of dual integrands is regular. It is important to monitor numerical stability when probing points close to such surfaces, as dual cancellations occur by cancelling large summands. We monitor this stability by testing the invariance of dual integrands under rotation of the spatial parts of the loop momenta integrated over. The more challenging loop integrals required a custom numerical stability rescue system that promotes the floating point arithmetic accuracy to quadruple precision when needed (which is about a factor 30 slower). We note however that the introduction of a complex contour deformation mitigates the numerical severity of dual cancellations.
Since we are mostly interested in verifying our method at this stage, we stress that no effort was made to finetune the integrator, sample statistics or loop momenta parametrisations. Sizeable improvements can be expected from considering techniques similar to the ones described in ref. [16]. Similarly to what was found in ref. [29], we observe that the Cuhre integrator offers significantly better convergence at one loop. However, we find that it performs much worse than Vegas at higher loops. For uniformity, we restricted ourselves to using the Vegas integrator for producing the results of table I. Our implementation is written in the Rust language, with Python bindings, and interfaces to the Cuba [37] library and Vegas3. 4 [38] for performing the adaptive Monte-Carlo integration.

IV. GENERAL KINEMATICS
The LTD expression of eq. 6, evaluated at external kinematics relevant for computing physical scattering amplitudes, typically features singular E-surfaces. The corresponding singularities then require a complex contour deformation of the spatial part of the loop variables, constructed so as to satisfy the LTD prescription associated to the surface, presented expanded at the first order in δ in (14). In sect. II A, we found that the sign of the imaginary part of the defining equation of E-surfaces reads: We now aim at constructing a contour deformation that satisfies the causality constraints implied by the iδ prescription. Given its parametrisation k C l = k l +i K l , l ∈ {1, ..., n}, one has that q j → q j + i κ j , ∀j ∈ b, and the imaginary part of every other propagator momentum can be expressed as a linear combination of { κ j } j∈b . This results in ∆ b,σ i acquiring an imaginary part, in the first order truncation of the expansion in j |s b ij | κ 2 j . For E-surfaces, this simplifies to which can be matched with (15) on individual Esurfaces by just setting One can observe that for every value of the loop variables where Thus, κ j ∝ − q j , ∀j ∈ b satisfies the prescription on arbitrarily many Esurfaces associated to the same loop momentum basis b, including on their intersection. Indeed, q j has positive projection on all non zero |s ij | v b i,j which might appear as summands in (17). This fails only if there exists an E-surface S b,σ l such that |s lj | v b l,j = 0, ∀j ∈ b, which would correspond to a pinched surface necessitating a soft/collinear regulator or subtraction. Finally, the intersections of two surfaces S b,σ a and Sb ,σ b with b =b would lie on dual cancelling surfaces, and are thus not singular in the sum of dual integrands. We stress that the above does not provide a complete recipe for building an overall continuous deformation direction K l satisfying all causal constraints and common to all dual integrands so as to preserve dual cancellations. This requires a (numerically efficient) strategy for interpolating between the deformation directions identified in eq. (18) for each group of E-surfaces E b j . Additionally, special care must be taken when setting the normalisation of the resulting deformation vector K l .
We conclude this section by presenting a first twoloop numerical result from applying LTD to a doublebox topology that requires a deformation around its 13 distinct E-surfaces. We set the external kinematics identical to those of the benchmark point chosen in ref. [9] and also report them in the ancillary material. Using Vegas3.4 with 105M Monte-Carlo samples, we obtained −5.877(55)·10 −14 , which stands within 1% of the analytical result −5.8973 · 10 −14 . We will provide a general and numerically efficient contour deformation for multi-loop LTD in an upcoming publication.

V. CONCLUSION
We derived a novel expression for multi-loop LTD, that involves taking as many on-shell conditions as there are loops. We demonstrated its potential for numerical integration by applying it to eight finite scalar multi-loop topologies. Additionally, we gave a first result of a contour deformation at two loops, showing that LTD can be used for computing integrals with physical kinematics as well.
Multi-loop LTD is a promising approach from both an analytic and a numerical perspective. One challenging possibility is to directly combine virtual and realemission unresolved degrees of freedom, as already explored at one loop in ref. [32], allowing for their joint integration with fewer or no counterterms. For numerical integration, LTD gives the advantage of reducing the number of dimensions to three per loop momentum, and confines singular surfaces one must deform around to a bounded region.
Our future work concerns extending the application of LTD to diagrams and loop amplitudes featuring (1) complicated overlaps of E-surfaces requiring a general contour deformation, and (2) UV and IR divergences, by designing local subtraction counterterms that lever-age known factorisation properties, such as the ones introduced in ref. [22]. [38] G. P. Lepage, J. Comput. Phys. 27, 192 (1978).  This section provides an example of the application of our LTD expression for a three-loop topology and it illustrates how the cut structure output by the Python script given in ancillary material is intended to be used.
Two loop momentum bases share the same cut structures if their set of signatures is identical, in other words if their defining loop momenta only differ by constant shifts. This implies that propagators sharing the same signature can be combined into a single list that we refer to here as a loop line, as done in ref [12]. A graphical equivalent of this procedure corresponds to constructing a reduced version of a loop diagram, where propagators with identical signatures are merged into a single loop line (see fig. 2(b)).
For a given reduced diagram (encoded as a list of loop momentum signatures), we can write the cut structure corresponding to each specific momentum basis (or equivalently: spanning tree) b i as a list Σ i containing as many elements as there are loop lines. At n loops, n of these elements are the members of σ bi and the remaining ones are set to 0 and serve to identify this particular spanning tree of the reduced diagram.
The cut structures of the diagram in fig. 2(a) (with 30 spanning trees) can therefore be obtained from the cut structures of the (simpler) reduced diagram in fig. 2(b) (with 8 spanning trees). As we see in fig. 2(c), both graphs have the same set of signature vectors.
For the reduced diagram in fig. 2(b) and when opting to close the contour of each energy integration in the lower-half complex plane, the Python script provided will return the following cut structures Σ i , each specifying all at once the spanning tree to consider and the cut energy signs corresponding to it: Σ 1 = ( 1, 1, 1, 0, 0), Σ 2 = ( 1, −1, 0, 0, −1), Σ 3 = ( 1, 0, 1, −1, 0), Σ 4 = ( 1, 0, 1, 0, −1), Σ 5 = ( 1, 0, 0, 1, −1), Σ 6 = ( 0, 1, 1, 1, 0), Σ 7 = ( 0, 1, 1, 0, 1), Σ 8 = ( 0, 1, 0, 1, −1). (A1) There are as many cut structures as there are spanning trees. For each cut structure, the element σ ∈ {±1, 0} at position i denotes that the loop line corresponding to s i is either cut with a positive/negative (σ = +1/ − 1) energy solution or not cut at all (σ = 0). When a loop line propagator is cut, one must consider one residue per propagator building this loop line. The energy solutions of these residues are then all equal up to constant shifts dictated by the energy component of the external momenta inserted along this loop line. When multiple loop lines are cut, all possible combinations of propagator cuts of each loop line must be considered. For example in the reduced diagram of fig. 2(b), cutting the loop lines q 1 , q 4 and q 5 (corresponding to the cut structure Σ 5 ) would yield 2 × 2 × 2 = 8 residues, while Σ 1 would yield only 2 × 1 × 1 = 2 residues. Alternatively, one can also opt to generate the cut structures of the original diagram of fig. 2(a) by treating each propagator independently. This would yield 30 cut structures, each in the form of a list containing 8 elements (one per loop line), each this time corresponding to a single residue.
We now make explicit the notation {q 0 j = σ b j E j } j∈b by writing, in the case of massless propagators, what are two of the four energy solutions defining the residues corresponding to the cut structure Σ 2 . These are the loop momentum energy configurations solving eq. (4) of the main text. A first residue stems from the energy solutions arising from cutting the propagators with momenta q 1 , q 2 and q 5 while a second one, for the same cut structure Σ 2 , comes from cutting the propagators with momenta q 1 + p, q 2 and q 5 and similarly for the remaining two residues.