Mathematical properties of nested residues and their application to multi-loop scattering amplitudes

The computation of multi-loop multi-leg scattering amplitudes plays a key role to improve the precision of theoretical predictions for particle physics at high-energy colliders. In this work, we focus on the mathematical properties of the novel integrand-level representation of Feynman integrals, which is based on the Loop-Tree Duality (LTD). We explore the behaviour of the multi-loop iterated residues and explicitly show, by developing a general formal proof for the first time, that contributions associated to displaced poles are cancelled out. The remaining residues, called nested residues as originally introduced in Ref. \cite{Verdugo:2020kzh}, encode the relevant physical information and are naturally mapped onto physical configurations associated to nondisjoint on-shell states. By going further on the mathematical structure of the nested residues, we prove that unphysical singularities vanish, and show how the final expressions can be written by using only causal denominators. In this way, we provide a mathematical proof for the all-loop formulae presented in Ref. \cite{Aguilera-Verdugo:2020kzc}.


Introduction
Quantum Field Theories (QFT) have shown to be one of the most successful theoretical constructions to describe the behaviour of Nature at sub-atomic scales. In order to extract reliable predictions from them, it is necessary to develop powerful methods inspired by mathematical ideas. The perturbative approach is one of these techniques, and nowadays it stands as the most important framework for high-energy particle physics. In this context, higher-order contributions in the perturbative expansion imply more precise and accurate predictions, with a reduced dependence on the factorization and renormalization scales. Then, it is crucial to compute such corrections in order to explore any small discrepancy with the highly-precise experimental data provided by high-energy colliders.
The computation of higher-order contributions in perturbative QFT is not straightforward. The main bottleneck is related to the calculation of virtual contributions, which involves dealing with multi-loop multi-leg Feynman integrals. In the context of current phenomenological studies [3][4][5][6][7][8], the calculation of scattering amplitudes involving oneloop Feynman integrals [9,10] is automated in different frameworks [11][12][13][14][15][16]. This great achievement is known as the Next-to-Leading order revolution. The same automation is not currently attainable for the evaluation of two-and, hence, multi-loop scattering amplitudes. The main obstacles, within standard and traditional approaches, rely on the reduction to the so-called master integrals, by Integration-by-parts identities [17,18], and the evaluation of multi-loop Feynman integrals (whose closed formulae are not all known as in the one-loop case).
In this paper, we mathematically elaborate on a promising alternative representation of multi-loop amplitudes and Feynman integrals. This approach is based on the Loop-Tree Duality (LTD) theorem [51][52][53], whose main aim is to open loop amplitudes into nondisjoint tree-level amplitudes. The genuine LTD theorem is valid in an arbitrary coordinate system. In specific applications, though, it is naturally defined in the Euclidean space of the spacial components of the loop momenta which is implemented by integrating out the energy components. Several calculations of scattering amplitudes at one- [54][55][56][57][58][59] and twoloop [60] have been provided within this formalism, as well as the numerical evaluation of multi-loop Feynman integrals up to four loops [2,61], which are based on the multi-loop LTD representation recently proposed in Ref. [1]. In this context, the LTD approach is currently drawing the attention as a novel tool aimed at overcoming many of the current bottlenecks, and alternative representations have been presented by other authors [62][63][64][65][66].
Besides the explicit calculation of loop integrals, LTD is the fundamental component of the Four-Dimensional Unsubtraction (FDU) framework [67][68][69], which aims at a simultaneous computation of real and virtual contributions directly in the four physical spacetime dimensions. Since loop integrals are re-expressed as phase-space ones, the infrared and threshold singularities are clearly identified in a compact region of the integration domain [70,71]. Thus, offering the possibility to implement alternative regularization strategies based on the local cancellation of physical singularities present in loop and realemission contributions [72][73][74][75][76][77][78][79][80].
In this paper, we will focus on the algorithmic methodology to derive compact dual integrand-level representations of scattering amplitudes. Following the ideas presented in a Ref. [1], we elucidate in more details all the mathematical concepts that are behind the multi-loop LTD representation. Furthermore, the conjectures originally presented are proven in the present work. In order to do so, we follow the application of the Cauchy residue theorem in succession at multi-loop level together with a short-hand notation that makes clear simplifications and proofs within this new framework. In this manner, this work represents the necessary mathematical basis to keep developing the program on the application of LTD to three [1,2], four loops [61] and beyond.
The outline of the paper is the following. In order to properly handle the ideas of the LTD theorem, we present in Sec. 2 an overview of the mathematical properties required to define the dual integrand for any scattering amplitude. For this purpose, we explain the concept of the iterated residue on a set of primitive variables, and we give some relations among them. In Sec. 3, we establish a useful notation for the multi-variable iterated residue of a meromorphic function, in order to simplify the symbolic manipulation of the expressions and take advantage of the properties behind this methodology. Once the notation is established, we perform the connection with the usual QFT formalism in Sec. 4. There, we motivate a practical definition of the topological classification of diagrams; in particular, we recall the multi-loop configurations presented in Refs. [1,2], Maximal Loop Topology (MLT), Next-to-Maximal Loop Topology (NMLT) and Next-to-Next-to-Maximal Loop Topologies (N 2 MLT). Later, we give some examples of the usage of this notation and computational algorithms to derive a formal proof of the all-order formulae presented in Ref. [1]. We put special emphasis on highlighting the reduction of complex topologies (i.e. N 2 MLT and higher) into nested convolutions of MLT ones, which also allows to explain the causal structure of the final compact results. Finally, in Sec. 6, we summarize this work and present an outlook of current developments, and potential future applications. Some detailed derivations are provided in the Appendices.

Multi-iterated residues and the Loop-Tree Duality
In this section, we establish the mathematical basis of the iterated residue approach for the multi-loop dual representation in the context of the Loop-Tree Duality formalism. Then, to begin the discussion, we will consider variables, functions and their pole structure, trying to identify the properties of the multi-variable residue. For this reason, we also introduce the physical concepts from the mathematical formalism, in order to appreciate the generality of the presentation.
We start by studying a multi-variable rational function f : R L → C. We will restrict the analysis to those rational functions involving quadratic polynomials in the denominator. Hence, we consider the following function where with β ∈ {−1, 0, 1}, k l ∈ R, x = (x 1 , . . . , x L ) ∈ R L is the primitive set of variables and y i ∈ C, for every i ∈ {1, . . . , m}, are the roots of the different polynomials that factorize the denominator. This functional form is sufficient to capture the essential mathematical structure of multi-loop multi-leg scattering amplitudes involving quadratic propagators and standard Feynman rules for the interaction vertices. As explained in Sec. 4, external momenta are encoded as shifts in the pole structure along the real axis. Keeping in mind the physical motivation, we restrict to the case γ i ∈ N which corresponds to having elements that strictly belongs to the denominator of the whole rational function. By convention, we take m ≥ L + 1, and with k L+1 ∈ R, an arbitrary number. Regarding the structure of the numerator, N , we will consider polynomial functions in R[ x] without any other restriction. It is worth appreciating that most of the proofs presented throughout this paper require very weak constraints on N , because this guarantees the generality of the approach and allows to extend their validity to almost any QFT.
Once the function f is described, we now move to the computation of its integral over the whole domain, R L . Explicitly, where the so-called primitive variables play the role of integration variables. Due to the rational structure indicated in Eq. (2.1), the most natural strategy to perform the integration consists in the application of Cauchy's Residue Theorem (CRT). In order to use this theorem, we must assume that f fulfils integrability on R L ; in other words, the integral in Eq. (2.4) must exist. Once integrability is guaranteed, we can: • use Fubini's theorem to change the integration order freely; • change integration variables (through linear combinations, re-scaling, etc.); • iterate the application of CRT.
At this point, the motivation for inquiring into this computation is related to the definition of multi-loop scattering amplitudes. As mentioned before, we can consider f as the integrand originated in any one-dimensional perturbative QFT and the integral in Eq. (2.4) as the associated amplitude. Moreover, a clever redefinition of variables will be enough to extend the treatment to QFT in an arbitrary number of space-time dimensions. In order to make use of CRT, let us take a look at the singularities of f . From the meromorphic structure of the function f (inspired on the integrands obtained in Feynman representation), it is seen that it only contains poles which arise as the solution of the equation, for the primitive variables. Notice that some solutions might have multiplicity higher than one, leading to multiple poles. The solutions of Eq. (2.5) for x 1 are given by, (2.6) for some l ∈ {L + 2, . . . , m}. Furthermore, when applying CRT, the integration variables are extended to the complex plane. Then, we can extend all of them simultaneously, or one-by-one. Also, we must take special care with those parameters entering in Eq. (2.4) that might displace the pole position in the complex plane. Thus, with the purpose of removing these ambiguities, we consider that: • the primitive variables are extended to the complex plane successively, but not simultaneously, which implies that when computing the residue for x i , we promote x i ∈ C but we keep x j ∈ R for any i = j; • every y i has positive real part and an infinitesimally small negative imaginary part.
Regarding the last point, we assume, with the purpose of defining the complex prescription. Again, this is inspired by physical concepts, and coincides with the customary +ı0 Feynman prescription introduced in QFT calculations. For the sake of simplicity, and to avoid overloading the notation, here and in the following, we drop the tilde, i.e. y i ≡ỹ i . More in details, the algorithmic procedure of the iterated residue begins with the promotion of the primitive variable x i to C, through the natural inclusion mapping i : Then, the residue of the latter function is computed along a contour included in the half-plane with Im(x i ) < 0 with the function Res : C C×R L−1 → C R L−1 . It is mandatory to say that the residue is well defined because the arguments of the function contain only one complex variable, x i . In this way, the iterated residue algorithm can be understood as an iterated application of the functor, and, thus, the iterated residue is represented by In the following, we will explore the consequences of these prescriptions. In particular, we will prove that some contributions vanish because of non-trivial cancellations related with the quadratic dependence of the denominators. As a consequence, we will explicitly show that the integral in Eq. (2.4) can be computed by just looking at the residues of specific poles.

Cancellation of residues from displaced poles
Now, we can choose one by one the integration variables and apply CRT on Eq. (2.4). We follow the natural order, i.e., x 1 first, x 2 then, and so on. The final result for the integral I in Eq. (2.4) is independent of the ordering, but intermediate expressions and integrand level results exhibit a non-trivial dependence on it. Then, to illustrate this behavior, we start by applying CRT in x 1 : we promote x 1 ∈ R → C and close the integration contour from the lower part of the complex plane, i.e.
where we assume that the function f fulfils the good-convergence hypothesis that allows integrability at infinity. The Heaviside theta function selects only those poles with negative imaginary part. Due to the fact that x i ∈ R ∀i = 1 and Im(y l ) < 0, it turns out that the only poles from (2.6) that contribute are for l ∈ {L + 1, . . . , m}. Then, Eq. (2.10) reduces to without any loss of generality. The next step consists in performing the second iterated integral in the primitive variable x 2 . This is the crucial step, because here it is necessary to recalculate the position of the poles in x 2 . After considering the poles in the variable x 1 , the position of some poles in the variable x 2 are displaced. At this point, and without giving the explicit formulae of the resulting integrand in each application of CRT, our purpose is to show that only the residues of specific poles contribute to the final result. We shall now select all the poles in the variable x 2 with negative imaginary part. Since {y i } i=1,...,m are the only parameters with non-vanishing imaginary part, the selection of the poles is done according to the specific combinations of y i 's that might appear as arguments of Heaviside theta functions.
In order to clarify the meaning of the last sentence, we present an explicit example. We assume that f takes the form and we compute the poles in x 1 . By selecting only those with negative imaginary part, we have 14) and the sum of residues is given by .

(2.15)
Then, by applying CRT on x 2 , the possible poles will be 16) and, from this set, we must only retain those with negative imaginary part. Performing the complete computation, we obtain where we can identify two kinds of contributions: • two terms associated to x 2 = y 2 and , (2.18) and Res( Res(f,

19)
• and one contribution which contains a non-trivial theta function, i.e.
[Res( Res(f, It is crucial to appreciate that, after the explicit computation, the sum of the iterated residues in Eq. (2.20) vanishes. These are the contributions associated to what we call displaced poles, whose position in the complex plane depend on the evaluation of residues in the previous variable. Moreover, in this example, the indices 1 and 2 for the integrated primitive variables were arbitrary; this implies that the result can be extended by induction to any number of primitive variables. Besides that, we would like to highlight that this property holds for any function f as described in Eq. (2.1), because it is a consequence of the quadratic pole structure. The master formula responsible of the cancellation is where a i and a ij are linear combinations of k l 's or y l 's, excluding the primitive variables x i and x j , and P is any meromorphic function without poles for the variables {x i , x j } in the location indicated in Eq. (2.21). Notice that Eq. (2.21) is valid for poles of arbitrary order. The proof goes through a direct computation of the iterated residue from the Laurent series of the function F (x i , x j ). A formal proof of this expression is available in Appendix A. This is the first mathematical proof of the cancellation of the displaced poles contributions and is the main result of this section. After each iteration the resulting function will have the form given in Eq. (2.22), so we conclude that Eq. (2.21) implies the cancellation of all the non-trivial Heaviside theta functions. Since this result holds step-by-step in the calculation, the final expression will also be free of residues associated to displaced poles. To conclude this example, we notice that when evaluating Eqs. (2.18) and (2.19) for the case L = 2, in Eq. (2.3), we obtain extremely compact expressions, namely, Moreover, we can appreciate that the denominator of the last expression in Eq. (2.23) only involve sums of y i 's. This observation is very important, since it restricts the presence of certain kind of singularities at integrand level that can be interpreted in terms of causality [2]. We will provide explicit all-loop order proofs that support the achievement of this causal structure in a very general family of topologies.
As the expressions obtained after the computation of the iterated residue are free of contributions from displaced poles, it is fair to ignore them. Thus, the remaining contributions are associated to what we call nested residues [2].

Towards a geometrical and physical interpretation
In the previous example, we observed two interesting properties, namely: 1. Cancellation of the residues from displaced poles: residues of poles whose imaginary part depends on different sign combinations of y i 's cancel.
2. Cancellation of non-causal contributions: only denominators involving sums of y i 's survive after adding up all the terms produced by the nested residues.
Regarding the first item, it means that terms proportional to θ (Im (y j − y i )) (2.24) cancel as successive evaluation of the corresponding residues leads to terms with opposite signs. These displaced poles are located in the upper or in the lower part of the complex plane depending on the specific value of the y i 's. On the contrary, the remaining contributions are those involving same-sign-combinations of y i 's. These poles always remain on one side of the real axis.
A geometrical interpretation of this cancellation is as follows. Let us start with the function . (2.25) The poles of the function in Eq. (2.25), in the variable x 1 , are located at as it is shown in Fig. 1. It is important to notice that x 1 has been extended to C, while x 2 is still considered as a real parameter. Thus the poles ±y 3 − x 2 are located along horizontal lines, depending on the value of Computing the residue of this function, closing the contour on the lower half-plane, the enclosed poles are each of which gives the corresponding residue, , .

(2.28)
Thus, it is obtained . (2.29) The poles of the first term are located at 30) and the poles of the second term are at Diagrammatically, the poles structure of each term in Eq. (2.29) are depicted in Fig. 2, where particularly the pole y 3 − y 2 is located somewhere inside the grey circle. This is because the imaginary part can be positive or negative, depending on the explicit values of y 3 and y 2 .
Since in this example we are dealing with simple poles, the computation of the residue is straightforward, although the overall interpretation that follows is also valid for multiple poles. In the first line of Eq. (2.28), the function must be evaluated in x 1 → y 1 and, in the second line, the evaluation is performed in x 1 → y 3 − x 2 − k 3 . On top of that, we notice that the factor in the denominator is symmetric under the transformation x 2 → −x 2 . As it is shown in Fig. 2, the location of the poles, except in the case y 3 − y 2 − k 3 , is symmetric with respect to the origin. This can be interpreted as a connection between the first and the second term through the transformation x 2 → −x 2 . This situation admits a graphical interpretation: performing the integration through the contour selected in the second term (on the lower half-plane) is equivalent to choose the integration contour on the upper halfplane for the first term. The subtlety here is that the pole y 3 − y 2 − k 3 will lead to a vanishing contribution because it appears inside both contours and thus, produces two contributions with opposite signs. To be more explicit, the imaginary part of the displaced poles could be positive or negative. On one hand, if the imaginary part of the displaced pole is positive, then it does not belong to the interior of any of the integration contours. On the other hand, if the imaginary part of the displaced pole is negative, through the transformation x → −x (which can be interpreted as a rotation in π) a new pole will appear in x 2 → −y 3 + y 2 + k 3 , leading to a relative minus sign between the contributions of the displaced poles. From the physical point of view, the cancellation of displaced poles is a consequence of the causal structure of multi-loop Feynman integrals, since these contributions do not have a representation in terms of cut diagrams. As discussed in previous studies [1,51,71], evaluating the residue of the integrand expression for a given loop diagram is equivalent to set on shell certain internal lines. The LTD representation is in fact independent of the loop momentum labelling and the ordering in the computation of the iterated residues. However, individual terms, associated to displaced poles, exhibit an explicit dependence on the loop labelling and ordering through the argument of non-trivial Heaviside theta functions. In consequence, they do not contribute to the final result, as we proved in the Sec. 2.1.
Up to now, we justified the cancellation of the residues of displaced poles relating them to unphysical contributions. In the same spirit, we can think about same-sign-combinations of y i 's as aligned contributions. The sign of the imaginary part of the poles is directly related with the energy flow of the internal propagators that are being set on shell. Thus, this means that only those contributions associated to a properly-aligned energy flow will remain. In the physics language, these are causal configurations and are directly related to the threshold singularities of scattering amplitudes in the context of the LTD approach, as discussed in Ref. [71].

Symbolic treatment of iterated residues
Once the mathematical basis of the nested multi-residue strategy was explained, we aim at simplifying the symbolic treatment of expressions. For this purpose, we deepen into the development of a physically-inspired notation that captures the main features of the iterated residue method and allows for a straightforward implementation.
We start from the conventions introduced in Sec. 2. When studying the analytic properties of the function given in Eq. (2.1), we can drop the information related to the specific variables and keep only the associated indices. In this way, we introduce the following notation: It is important to notice that there are no common zeroes between F i (i) −1 and F j (j) −1 when i = j. Thus, by using the definitions of Eq. (3.1), the original function in Eq. (2.1) can be rewritten as From now on we consider the case N = 1, but we remark that most of the results remain unchanged as the concept of the residue does not depend on the specific structure of the given function. Instead, it depends only on the pole structure (i.e. in the zeroes of denominators). This notation has some interesting properties which allows to simplify the presentation of explicit results. In particular: 1. Having two or more indices together in the k-th argument is interpreted as summing the variables associated to those indices inside the propagator corresponding to y k . For instance, 2. Having a bar above a given index means that the associated variable is inverted, 3. Having sub-indices within a given argument means that associated y-parameters are added or subtracted according to the bar convention. For example, 4. Having 0 as one of the arguments represents that the corresponding primitive variables where replaced by the combination of y-parameters present in the sub-indices. Also, if 0 is in the k-th argument, and if it has as sub-index (k), it means that it has been computed the residue in the poles within the associated factor. E.g., which corresponds to computing the residue in x 1 = y 1 . It is important to emphasize that, for a function such as F (1, 2, 12), the sub-index 3, representing a term y 3 , is different from the sub-indices 12, representing the sum y 1 + y 2 .
Some immediate properties of the quadratic structure of the functions F i are that F i (j k ) = F i (j k ) and F i (j k ) = F i (j k ), and thus, it is possible to fix the overall sign of a given variable. Also, it is straightforward that k = k and kk = 0. This two properties allows us to compute efficiently the nested residues.

Efficient residue computation
At this point, we are interested in computing efficiently the residues of a generic function with the form F (1, 2, . . . , m). Moreover, we aim to use algebraic and symbolic properties to avoid computing unnecessary terms which cancel in each iteration of the residue. For illustrative reasons, we focus the discussion on the case of functions with the particular form, In the following, we will use the equivalent notations L with the purpose of shorten the presentation of results, wherever the expressions are unambiguously defined. The poles with negative imaginary part of F , in the complex variable x 1 , are Whence, for the computation of the residue at x 1 = y 1 , we obtain, whilst the residue in the other pole is given by In general, after computing a residue, we obtain different functions of the form F i (i), F i (j), F i (j k ) and F i (j k ). The set of negative imaginary part poles associated with these functions are As already noticed in Sec. 2, the iterated residue might lead to expressions whose poles are not always within the integration contour (i.e. with negative imaginary part); the so-called displaced poles. For instance, y k − y i + k L+1 might have a positive or negative imaginary part, depending on the specific values of the y-parameters. Thus, it is necessary to impose the condition Im(y k − y i ) < 0 when computing the residue with the function θ(ik) := θ(Im(y k − y i )). (3.12) Hence, Remarkably, from Eq. (3.13), we directly appreciate the cancellation of the contributions associated to displaced poles, due to the appearance of a relative minus sign. In consequence, this cancellation becomes explicit for the second and subsequent iteration of residues. The formal proof presented in Appendix A is based on identifying the contributions with opposite signs that cancel among them. Notice that similar cancellations have been observed in [64,65] by considering poles with positive and negative imaginary parts for specific configurations without a formal general proof as presented in this paper.
In order to clarify the notation, we present an explicit example. For instance, where the arrow represents the computation of the residue of function on the left. In the first line, the residue in the variable x 1 originates two terms. Then, in the second line, we identified all the poles in x 2 associated to the expression in line 1, and we computed the residues. By putting the sub-indices and subtracting from the main index, we identify the variable in which we apply CRT and the pole where we evaluate. The extra contribution F (0 (1) , 0 13 , 0 (3) ) has not been considered because it corresponds to a displaced pole. A practical way to see this procedure is as follows. The computation of the first residue in Eq. (3.14) corresponds to the poles x 1 = y 1 and x 1 + x 2 = y 3 + k 3 , or, in this notation, corresponds to 1 = 0 1 and 12 = 0 3 . These two index equations are equivalent to 1 1 = 0 and 12 3 = 0, respectively, and, as these expressions are 0, they can be added or subtracted in other arguments in the corresponding iteration of the residue. Thus, for the computation of the first residue, in the case of the pole x 1 = y 1 , it can be obtained 12 = 12(1 1 ) = (11)2 1 = 2 1 . Analogously, for the pole Finally, as 2 3 becomes the first argument of the function, it represents a factor of the form It is important to highlight that, through the computation of the iterated residue, the displaced poles can clearly be identified in two cases. In the first case, the displaced poles can be seen as those arriving from arguments with at least one sub-index without bar. For instance, in Eq. (3.14), after the computation of the first iteration of the iterated residue, with respect to the variable x 1 , in the first term, F (0 (1) , 2, 2 1 ), the poles associated to the third argument 2 1 are one positive-imaginary-part pole (which is outside the integration contour) and one displaced pole. The second case corresponds to the arguments with all sub-indices with bar. In this case, there is one negative-imaginary-part pole and one displaced pole. Then, the displaced pole is identified as the one leaving the argument with at least one sub-index without bar and at leas one sub-index with bar. For instance, in Eq. (3.14), after the computation of the first iteration of the iterated residue, with respect to x 1 , the second term F (2 3 , 2, 0 (3) ) has its first argument 2 3 , and thus contributes with one displaced pole (located in x 2 = y 3 − y 1 + k 3 , or, in this notation, 2 = 0 13 ) and a negative-imaginary-part pole (located in x 2 = y 1 + y 3 + k 3 , or, in this notation 2 = 0 13 ). In this manner, the nested residue can be computed directly by considering just the poles associated with the arguments without sub-indices and the poles associated with the arguments with all its sub-indices with bar, whenever the index of the integration variable does not have a bar.

Recursive representation with nested residues
In Refs. [1,2], we showed very compact formulae for some Feynman diagram topologies at all-loop orders. We provide in this paper formal proofs of the validity of these expressions, by taking advantage of the notation previously introduced with the purpose of unveiling the recursive relations that naturally manifest when computing the nested residues. This will lead to inductive proofs of the beforehand mentioned formulae.
In the following, we explore some relations among nested residues to identify the potential recursive structures. We sequentially calculate the residue in the primitive variables, and simplify the result of each step to find the dependence on the number of iterations. Thus, we infer the functional forms that we obtain after the i-th iteration and proof the inductive step.
So, let us study some examples to find the recursions. The simplest case corresponds to the application of the iterated residue to a function F (1, . . . , L), with independent arguments. This is the case of which corresponds to the class of factorizable functions and does not deserve further comments. Here, we discuss about the non-trivial case of non-factorizable functions, and the most symmetric example is found when there are only two dependent arguments in each step of the iterated residue, as for the function F (1, . . . , L + 1). We start by noticing an interesting property when computing the first i-th nested residue for this function, where F (1, . . . , L − i) was factorized because it does not depend on the primitive variables that are involved in the computation of the iterated residues in the last i-th variables.
Whence, it is straightforward that, after the computation of all the iterated residues, we obtain . . , 0 (L+1) ). (3.17) In Appendix B, we provide a formal proof of these relations for simple poles, by using a more physically-inspired notation and, in Appendix D, we prove the generalization for multiple poles. Again, we recall that the arrow indicates that the expression in the r.h.s. corresponds to the nested residue of the expression in the l.h.s. Another important relation can be found for functions of the form F (1, . . . , L + 2), where we defined L + 2 ≡ 12 = −1 − 2. In this case, after a direct computation of all the iterated residues and reordering the result, we obtain The nested residues lead to two terms, with a strong resemblance to the factorization formulae presented in Ref. [71]. The different contributions are characterized according to the poles considered for the residue computation. In the following, we will explain better this separation, although we defer to Sec. 5 the physical interpretation in terms of on-shell internal propagators.
Let us extend the results for more general functions. For a given function F (1, . . . , m), if {1, . . . , ρ} are the indices of the primitive variables appearing in three or more arguments then, whenever γ i = 1 and the k-parameters vanish for i > ρ, direct computation of the first L − ρ residues of the function F leads to where (3.20) The right hand side of Eq. (3.19) encodes the singular behaviour of the function in the left hand side and the function defined in Eq. (3.20) represents an auxiliary propagator which summarize the information associated to the sets ρ + 1 through L + 1. From a physical perspective, this expression plays the role of a modified propagator with an alternative on-shell condition. Also, it can be thought as a consequence of applying momentum conservation; we will return to this point later, in Sec. 5. Notice that it is enough to show the validity of these expressions for simple poles (the formalization of this claim is given in Appendix D). The generalization to the non-vanishing k-parameters case presents no extra difficulty and is delayed to Sec. 5.
Thus, the r.h.s. of Eq. (3.18) can be expressed as where After a direct computation of the iterated residue, we end up with where it is understood that

(3.25)
In the first case, it is important to notice that it is not possible to consider simultaneously the pole associated to the function F L+2 , so that the function F L+2 (L + 2) = F (L + 2) is factorized.
We would like to highlight that the great simplification from Eq. (3.18) to Eq. (3.23) is due to the fact that all the information of the function F (3, . . . , L + 1) is encoded within (3.28) With this identification between the (L − ρ)-th iterated residue of the function F (ρ + 1, . . . , L+1) with multiple-poles or non-vanishing a-parameters, and the function F (ρ+1)...(L+1) which is its simplification for simple poles and vanishing a-parameters, it is possible to interpret the arguments {ρ + 1, . . . , L + 1} as a single function with only one factor in the denominator of the form (x 2 − y 2 ) with multiplicity 1.
All this discussion is needed to consider the case of a more complex function of the form F (1, . . . , L + 3), where we define L + 3 ≡ 23 = −2 − 3. In this case, we factorize the sets 4 through L + 1 and Eq. and  and for the rest of the arguments, it is the sum of functions of the form F (4, . . . , L + 1) where the iterated residues have been computed and one of the poles associated to the arguments {4, . . . , L + 1} has not been taken into account. The results obtained in Eqs. (3.23) and (3.30) are remarkable because they give a relation between a family of functions with certain complexity and convolutions of simpler ones. This association can be understood in a diagrammatic way, by using the concepts and ideas taken from QFT. We will discuss that in the following section, and apply these relations to provide recursive proofs of some physical results in Sec. 5.
We would like to highlight that the discussion of this section was independent of the numerator of the integrand f and independent of the order in which the residues of the poles of F are evaluated, which ensures that it is general enough to allow a straightforward application to scattering amplitudes computations.

Nested residues for scattering amplitudes
In the context of perturbative QFT, we are interested in computing scattering amplitudes, and, in particular, their higher-order representations. These representations are given in terms of loop Feynman diagrams, where internal virtual particles circulate as quantum fluctuations. A diagram with L loops posses L independent or primitive loop momenta, { i } i=1,...,L , that define the integration space. The momenta flowing through the internal lines can be grouped into sets, in such a way that all the momenta inside a set s are of the form q is = s + k is , with s and k is linear combinations of primitive loop and external momenta, respectively. Because of momentum conservation, the number of momentum sets, n, is always larger than the number of loops, n ≥ L + 1 for L ≥ 2 non-factorizable Feynman diagrams. 2 Likewise, the loop diagrams involve Feynman propagators, whose dependence on the loop momenta settles the pole structure of the whole amplitude. Thus, we introduce the scalar Feynman propagator, where q i represents the momentum flow through this line, m i is the mass of the particle and q (+) is the corresponding positive on-shell energy. The +ı0 prescription is crucial for establishing the location of the poles in the complex plane, and thus defining the physical modes for the on-shell states. In order to write down explicit representations for scattering amplitudes, we need to use the corresponding Feynman rules. In general, an L-loop amplitude with N external particles is given by where N is a function given by the Feynman rules of the theory that depends on the loop and external momenta, and is a product of Feynman propagators over the union of the n momenta sets, allowing arbitrary positive powers γ i for each line. As usual, the integration measure is defined as where {k m,0 } are linear combinations of the energies of external momenta and N is a polynomial in the loop energies. At this point, the connection with the notation introduced in Secs. 2 and 3 is straightforward: Eq. Regarding the short-hand notation introduced in Sec. 3, we identify a Feynman propagator associated to a line i ∈ s, G F (q is ), with F is (i s ). The nested residues correspond to the so-called dual amplitudes, as defined in Eqs. (5)-(6) of Ref. [1]. Explicitly, we establish the connection (4.7) We anticipate that this result justifies the so-called MLT formulae presented in Ref. [1], and shall be explained in more detail in Sec. 5. Finally, we would like to make a comment on Eq. (4.6). After the application of the iterated CRT, the original loop amplitude will involve only integrals in the spatial components of the loop momenta, i.e. i , which are inside the definition of q (+) i,0 . The integration space is now Euclidean, instead of the original Minkowskian one. This fact, together with the compact form of the dual representation, points towards a more efficient numerical implementation within this formalism, as we already tested in Ref. [2].

Topological families
Scattering amplitudes can be classified according to their internal momentum flow, which translates into specific topological structures for the associated Feynman diagrams. In Ref. [1], we introduced a systematic classification scheme of multi-loop topologies, which includes specific families of diagrams with arbitrary number of loops.
Given an L-loop diagram (L ≥ 2) with n ≥ L + 1 sets of internal propagators, we define the topological complexity as k = n − L. In this way, the Maximal Loop Topology (MLT), which is the most symmetric configuration, has topological complexity k = 1 and the Next-to-Maximal Loop Topology (NMLT) has topological complexity k = 2. In general, a N k−1 MLT diagram at L loops has topological complexity k, and we will denote it N k−1 MLT(L).
With all these definitions in mind, we proceed to present explicit results for MLT(L), NMLT(L) and N 2 MLT(L) configurations in the following sections, focusing on their recursive structure and the decomposition into convolutions of lower-complexity topologies.

Selected results for topological families
Multi-loop scattering amplitudes with an arbitrary number of external legs are objects that involve integrands with a structure that can be properly described in terms of the functions F i (j) defined in Sec. 3. In this section, we make use of their properties shown above, in order to highlight their recursive structure. The most symmetric loop configuration in any QFT can be encoded into the MLT(L) diagram which is given by and is graphically depicted in Fig. 3.  For the moment, we have not discussed anything regarding the structure of the numerator in Eq. (4.6), however, as discussed in Sec. 2, the treatment at integrand level is independent of the explicit structure of the latter. It is straightforward to notice that the symbolic handling of the expressions relies on the iterated application of CRT, that only requires to indicate the pole location, without making use of the explicit functional form of the numerators. Hence, for the sake of simplicity, we restrict the following discussion to the case N = 1, since all our dual representations can be straightforwardly generalised to any numerator. We can also restrict the demonstrations to vacuum diagrams, i.e. those without external particles, because they contain sufficient information regarding the loop-momenta dependence of each internal set of propagators. The generalization to loop configurations with an arbitrary number of externa particles is achieved by implicitly considering the sum over nested residues within each set of propagators.
In order to find the LTD realization of Eq. (5.1), we just need to interprete Eqs. (3.1) and (3.17) in terms of Feynman and dual propagators as, where we introduce the dual representation of an MLT(L) topology, and the arrow is used to indicate that the expression in the r.h.s. is the result of applying CRT to the original amplitude 3 . As pictorially depicted in Fig. 4, this formally proves the validity of the MLT(L) formulae presented in Ref. [1]. In the specific case of MLT(L) topologies with single powers and one propagator per loop set, we formally proof the formulae presented in Ref. [1]. After summing over all the dual L + 1 contributions, and applying the results given in Appendix C, Eq. (5.2) collapse to the extremely compact and causal expression This is a multi-loop generalization of Eq. (2.23). It is worth mentioning that Eq. (5.3) is a consequence of the algebraic properties of the nested residues, and the same strategy can be applied in order to show the explicit causal representations for more complex topologies as exhibited in Ref. [2].

NMLT(L) and N 2 MLT(L)
The MLT configuration is sufficient to describe any two-loop scattering amplitude, but new mathematical structures appear at higher orders. Starting at three loops, we also need to consider the NMLT and N 2 MLT topologies. In fact, the NMLT topology is described as a subtopology of N 2 MLT, which is the master topology at three loops. This constitutes a clear and powerful classification scheme towards an efficient computation of higher-order amplitudes, since these new topologies involve new loop momenta linear combinations. The NMLT(L) and N 2 MLT(L) are respectively characterised as follows; where we include two extra sets, i.e. {L + 2, L + 3}, with the aim of describing all the possible momenta configurations. This expression can be understood in terms of loop configurations of lower topological complexity as it is shown in Fig. 6. As it was anticipated in Sec. 3, the convolution symbol is not a pure factorization. We would like to emphasize that it implies the use of the on-shell conditions to express all the off-shell variables. We observe that the NMLT(L) displays two contributions in terms of MLT configurations. Explicitly, the first contribution is a convolution of an MLT(L − 2) with an MLT(2) diagram, and the second term is a convolution of an MLT(L − 2) diagram, all of its propagators are on shell and are reversed, with an MLT(1) one.
Finally, let us now draw our attention to the N 2 MLT(L) configurations. It is worth appreciating that this mathematical object contains the highest topological complexity at three-loop level. In order to describe it, we need to add the set L + 3 = 23, which is the set that can only contain combinations of the loop momenta 2 and 3 . The generated vacuum topology, or Mercedes-Benz like-diagram, is depicted in Fig. 7.  In order to find the dual expression of Eq. (5.5), we proceed to apply the CRT consecutively and, after computing the first residue and taking into account Eq. This decomposition is graphically described in Fig. 8, where the convolution symbols have the same interpretation as before. The brackets notation used for G D ( 1, L+3 , 2, L+2, 3 ) corresponds to the insertion of external momenta in specific internal lines. For instance, let us take a look at the second line of Fig. 6, where the dot between 1 and L + 3 represents the insertion of an external particle with momenta (L + 3)1. It is also important to remark that we obtain a dual expansion with two contributions where one of the terms is a convolution of an MLT(L − 3) configuration with an NMLT (3), and the other term consists of one MLT(L − 3), all of its internal momenta set on shell and reversed, and an MLT(2) configuration with two propagators in two internal lines. The decompositions of NMLT(L) and N 2 MLT(L) topologies shows explicit recursion relations involving configurations with lower topological complexity. Therefore, from the above studies and by an iterated application of the decomposition presented in Figs. 6 and 8, we can notice that any N k−1 MLT (for k ≤ 3) can be cast in terms of MLT configurations. Besides, it is interesting to point out that recent studies that include master topologies at four-loop also present the same behaviour [61]. Therefore, an extensive study of MLT configurations and their convolutions, as carried out in the present paper, is sufficient to understand the behaviour of any L-loop amplitude with any number of external legs.

Higher topological complexity and causality
The ideas presented in this paper can be generalized to scattering amplitudes with an arbitrary topological complexity. This is due to the fact that the algorithm for the computation of the nested residue does not depend on the number of loops nor the topological classification of the diagram. Moreover, the algorithmic procedure is the same whether or not external particles are present. Figure 9. Graphical factorization of a multi-loop topology with a MLT insertion. The gray blob represents the subtopology with specific topological complexity. The diagram on the r.h.s. of the arrow represents the opening of the MLT subtopology.
In order to deal with an amplitude with arbitrary topological complexity, it is useful to express Eq. (3.19) as it is shown in Fig. 9, where α = 1 . . . ρ * plays the role of the internal line L + 1 in the minimal Feynman diagram with the same topological complexity. In other words, the line α is equivalent to unifying all the MLT-like insertions of the original diagram In this figure, the topological complexity of the diagram has been isolated inside the blob J, and the remaining MLT-like part of the diagram is simplified with the results of this work. This is important because, after the computation of the iterated residue (equivalently, after a partial opening the diagram), the presence of external particles attached to the vertices isolating the topological complexity can be thought as merged into a new internal line, whose momenta flow is determined by momentum conservation. We would like to highlight that Eq. (3.19) nor Fig. 9 are final results, since still remains the nested residues with respect to momenta 1, 2, ..., ρ and α have to be computed.
Also, for interactions of external particles with the internal lines of the subdiagram J, the analytical structure remains untouched: all the information regarding the external lines is codified inside the factorized contribution associated with J. Thus, we can only think about a different J , clearly with additional poles, but the factorization formula remains the same. Some works have been developed for the study of the N 3 MLT(L) and N 4 MLT(L) [61]. There, it is seen that the computation of the nested residue for topological complexities 4 and 5 yields to representations analogous to Eqs. (5.6) and (5.7).
If external particles are attached to an internal line of the MLT subtopology, the sum over the nested residues of all the propagators that belong to the corresponding set is required. All the poles within each set have imaginary parts of the same sign. In this situation, the generalization of the results derived from Eq. (3.21) becomes straightforward: the nested residue takes the form of a sum of as many copies as propagators are in the same internal line, where the poles are shifted one to another by a real number. This configuration is depicted in Fig. 10. The application of the ideas presented in this paper will be useful when a realistic scattering process is considered; where it will become mandatory to study the consequences of having a polynomial in the energy components of the loop momenta q i,0 as numerator. For the purposes of this work, it is not necessary to make an explicit example to claim that when the numerator is not identically 1, the results presented along this document are still valid if the numerator is a meromorphic function in every energy variable. And, since this is the case for a Feynman integral (the integrand is always a rational function of the energy of the loop momenta), then these results stand for any QFT. In addition, this approach can be used to obtain the causal structure of an arbitrary topological class of diagrams, since the nested residue leads, in a natural way, to sums of on-shell energy of the internal particles, q (+) k,0 (similar to the expressions of Appendix C), avoiding the non-physical threshold singularities. These causal structures makes it easier to localize the physical thresholds, as they will play an explicit role within the causal denominators.
Finally, we would like to make a brief comment about causality and the location of physical thresholds. As discussed in Ref. [2], when adding up all the dual contributions, the resulting expression is written in terms of causal denominators, {λ ± i }. These variables represent sums of on-shell energies and combinations of the energy of the external particles. The number of causal denominators depends on the topological complexity and on the number of external particles. However, their functional form and explicit dependence on the number of loops can be inferred from the causal denominators present in the associated vacuum diagram. Thus, the proofs provided in this article allow to ensure the validity of the all-loop order formulae presented in Ref. [2].

Conclusions
The computation of scattering amplitudes at higher-orders in the perturbative expansion is a very challenging task, specially for multi-leg processes. Even though several highly innovative and groundbreaking techniques were developed in recent years, automation of multi-loop scattering amplitudes still remains a frontier problem. In this respect, the Loop-Tree Duality offers alternative representations of generic scattering amplitudes at integrand level, which have many potential advantages over the customary Feynman representation.
In this paper, we deepened into the mathematical aspects of the multi-loop construction presented in Ref. [1]. Firstly, we provided a rigorous definition of the multi-iterated residue computation using a generic test function and exploring the consequences of the prescription introduced. We found that it is not necessary to keep the residues of all the poles whilst performing the iteration, since some of them cancel. These displaced poles are associated to non-physical contributions, which cannot be interpreted in terms of cut diagrams. This allows to redefine integrand representations in terms of the so-called nested residues which are only those related to physical contributions.
Inspired by the mathematical properties of the nested residues, we defined a closed notation to achieve an efficient symbolic handling of intermediate expressions. Following this approach, carefully explained in Sec. 3, the residue computation is performed by lowering indices, which contain the relevant information about the pole location. In addition, the nested application of the Cauchy residue theorem leads to recursive structures, whose behaviour is accurately captured by this short-hand notation.
Once the formalism presented here was provided with a physical meaning, in Sec. 4, we managed to apply it to specific benchmark amplitudes. In particular, we introduced a formal definition for the topological complexity of families of loop diagrams. We used these concepts and the operational methodology presented in Sec. 3, to inquire into the recursive structure of MLT(L), NMLT(L) and N 2 MLT(L) topologies. In this way, we provided the ingredients required for an inductive proof of the all-loop order formulas presented in Refs. [1]. Moreover, we showed that the recursive nature of the computations leads to an explanation to the causal behaviour of the compact formulae found in previous papers [2]. As detailed in Sec. 5.2, the ideas that we developed can be straightforwardly applied to any multi-loop multi-leg amplitude, independently of their topological complexity [61].
In summary, we exhaustively focused on the analysis of the mathematical structures behind scattering amplitudes by applying the LTD framework. From the formal properties that we found, we obtained valuable information for proving explicit all-order formulas, and also to efficiently perform the symbolic handling of the expressions. This knowledge allows to reach higher-perturbative orders, thus opening an interesting path for more precise theoretical predictions.

A Cancellation of residues from displaced poles
The cancellation of the residues from displaced poles, defined in Sec.2, is guaranteed by the following: Lemma: Let P (x i , x j ) be a meromorphic function in both variables x i and x j whose poles are not located on Then, the iterated residue in each of the explicit poles satisfies Proof : If the shifts x i = x i − k i and x j = x j − k ij + k i are performed, the function F can be rewritten in the form Without loss of generality, this is also equivalent to consider k i = k ij = 0. The function in Eq. (A.3) has two explicit poles of order γ i and γ k within the half plane Im(z) < 0. Thus, the function F has an expansion of the form If the last factor of the right hand side of Eq. (A.4) is rewritten in the form x i + x j − y k = (x i − y i ) + (x j − y k + y i ), and if the sum over r k is split into negative and non-negative values, it is obtained To compute the first residue of the function in Eq. (A.4), for {x i , y i }, it is enough to take the coefficient of the term with the factor (x i − y i ) −1 in the expansion of Eq. (A.5). Afterwards, to obtain the second residue, for {x j , y k − y i }, we select the coefficient of the term with the factor (x j − y k + y i ) −1 . For the second sum, the second condition is never satisfied because 0 ≤ s ≤ r k and such a factor demands the condition s = r k + 1. For the first sum, the second condition is obtained for r k + s = 1, and as 1 ≤ r k ≤ γ k and 0 ≤ s < ∞, there is just one term satisfying this condition, with s = 0 and r k = 1. As s = 0, the first condition is satisfied for r i = −1. Hence If in the function in Eq. (A.4), we rewrite the second factor as x i − y i = (x i + x j − y k ) − (x j − y k + y i ), and if the sum over r i is split into negative and non-negative values, we obtain Again, the iterated residue of this expression is the coefficient of the terms proportional to (x i + x j − y k ) −1 and (x j − y k + y i ) −1 . For the first sum, this conditions are satisfied for r k + s = −1 and r i + s = 1. However, as 1 ≤ r i ≤ γ i and 0 ≤ s, the last condition is fulfilled only for r i = 1 and s = 0. Thus, the first condition is expressed as r k = −1. For the second sum, the first condition holds, but the second condition shall be expressed as r i − s = −1 so that s = r i + 1. However, it is given that 0 ≤ s ≤ r i and then this sum does not contribute to the residue. Thus, It is then concluded that If we then restore the original variables that are shifted by k i and k ij with respect to x i and x j , we arrive to the expression we wanted to demonstrate (A.10)

B Proof by induction of the multi-loop MLT(L) representations
This Appendix presents a formal proof of the dual representation of MLT(L) in terms of nested residues (Eq. (3.17)). The proof is given by induction on the number of computed residues through the iterated residues algorithm. Here, we start by analysing the dual representation of a scalar MLT(L) diagram with one propagator for each set. The original integrand in the Feynman representation is given by I In order to prove the cancellation of the contributions of the displaced poles in each iteration of the iterated residues by mathematical induction, we assume that the function obtained after computing the first i iterated residues (for the last i variables) is given by where it has been factorized out the Feynman propagator G F (1, . . . , L − i) as it depends on independent primitive variables. Then, the set of poles of the function given in Eq. (B.3) with respect to the variable q L−i,0 is given by while all other poles are displaced poles, we should select only the residues of the nondisplaced poles with negative imaginary part. Following with the next nested residue, we get, It is worth to say that the proof of Eq. (B.3) is general enough to cover the case involving an arbitrary topological complexity. This is because we can isolate the highertopology structure inside the factor G F (1, . . . , L − i − 1), and proceed as described. Thus, Eq. (B.3) can be applied to Feynman diagrams with higher topological complexity and any number of loops.

C Causal rearrangement of nested residues
Let {1, . . . , ρ} be the family of sets with loop momenta appearing in three or more sets, and let ρ + 1 ≤ i ≤ L. After the computation of the i-th iterated residue of the integrand of a general Feynman integral, It is important to notice that the function in Eq. (D.2) is not a physical propagator. This is the reason we call it auxiliary propagator. It is worth to notice that this function is biyective, so that the inverse image ψ −1 is a function.
After the computation of k-th iterated residues with respect to L, L − 1, . . . , L − k + 1 to every element of F L,m it is obtained the subset Res k [F L,m ] of C (R L−k ) . Finally, let us define the operator Using the identity mapping in N m , id : N m γ → γ ∈ N m , we show that the algebraic diagram in Fig. 11 commutes. In Ref. [2], this transformation is used to relate the causal structure obtained in Ref. [1] for the MLT(L) with simple Feynman propagators and the expression for a double pole in one of the internal sets. In this work we use it in order to generalize the application to an arbitrary topological complexity, putting on the surface the sufficiency of the simple poles case. For a general topological complexity diagram, it can be written Whence, the computation of the iterated residue of the integrands in both sides of Eq. (D.9) with respect to the same variables in the same order leads to the commutation of the algebraic diagram in Fig. 11.