Lotty – The loop-tree duality automation

Elaborating on the novel formulation of the loop-tree duality, we introduce the Mathematica package Lotty that automates the latter at multi-loop level. By studying the features of Lotty and recalling former studies, we discuss that the representation of any multi-loop amplitude can be brought in a form, at integrand level, that only displays physical information, which we refer to as the causal representation of multi-loop Feynman integrands. In order to elucidate the role of Lotty in this automation, we recall results obtained for the calculation of the dual representation of integrands up-to four loops. Likewise, within Lotty framework, we provide support to the all-loop causal representation recently conjectured by the same author. The numerical stability of the integrands generated by Lotty is studied in two-loop planar and non-planar topologies, where a numerical integration is performed and compared with known results.


Introduction
The predictions of physical observables at high energies through the perturbative approach of Quantum Field Theories has drawn the attention of the scientific community when trying to understand in more details the physics at experiments. In fact, a support of the theory on the experimental side needs to be provided in view of the synergy with the collider machines, such as LHC at CERN [1], and the future colliders that aim higher luminosity [2][3][4][5][6][7][8][9][10][11]. Thus, to calculate these observables, we rely on scattering amplitudes, which turn out to be the backbone in these predictions.
Despite the plethora of diagrams that need to be considered when going to higher orders, new approaches, based on clever mathematical ideas, are currently being proposed to increase the efficiency and reduce the computing time in the several calculations [39][40][41][42][43][44][45][46]. Then, with these ideas and methods at work, the clear motivation is extending the knowledge of NLO to higher orders, N...NLO. However, the transition has not been straightforward and many subtleties that were not present at NLO have appeared. For instance, the evaluation of multi-loop Feynman integrals is in fact an open problem, since, differently to the one-loop case, a complete basis of master integrals is not known and yet need to be computed [47][48][49][50][51]. Nonetheless, calculations, at two-loop level, including five external partons, have started to appear displaying relevant results for the physics at LHC [52][53][54][55][56][57][58].
Hence, to provide an alternative approach to deal with the numerical evaluation of multi-loop scattering amplitudes, we elaborate on the loop-tree duality (LTD) formalism, originally developed in Refs. [59][60][61][62][63], and successfully employed in applications at one [64][65][66][67] and two [68,69] loops. The main idea of LTD, together with its motivation, is to open loops into connected trees, in such way that all contributions, independently of the loop order, are equally treated. This is, for instance, the case of the four-dimensional unsubtracted (FDU) scheme [70][71][72], which treats virtual one-loop with real tree-level diagrams as phase space integrals. This is carried out, in order to locally cancel, before integrating, infrared (IR) and ultraviolet (UV) singularities. Apart from the interplay of LTD with FDU, other formalisms that aim at having local cancellations of IR and UV singularities as well as a purely four-dimensional representation of the space-time dimension are currently being extended to have a complete automation at N 2 LO [70][71][72][73][74][75][76][77][78][79]. For recent reviews, we refer the reader to Refs. [80][81][82].
In the spirit of harnessing from the knowledge of LTD, a novel formulation of the latter was recently provided in Ref. [83], where the approach to obtain a dual representation of multi-loop Feynman integrals, within the LTD framework, allowed for a complete automation and analysis regardless of the loop order, e.g. the works of Refs. [84][85][86][87] perform an analysis of up-to four loop topologies. Besides the generation of dual integrands, we have recently observed that LTD satisfies a representation free of unphysical singularities when all dual integrands are summed up. For the latter we refer to the causal representation of Feynman integrands, since these integrands are expressed in terms of causal (physical) thresholds. Explicit results, together with numerical evaluations, for topologies up-to three [84] and four [85] loops were provided by reconstructing their analytic expression over finite fields [88][89][90][91]. To this end, we profited on the available tools and made use of the FiniteFlow framework [90]. We remark that, apart from our novel representation, alternative approaches to the LTD formalism have been addressed in Refs. [92][93][94][95][96][97][98]. In particular, let us remark the implementation of Ref. [97] with its stand-alone code cLTD.py.
In view of the above-mentioned features of the novel formulation of LTD, it is desirable to provide a pedestrian code to study in more details how the dual integrands and, therefore, the causal representation is generated within this framework. In fact, the algorithm proposed in Ref. [83], considered firstly as the iterative application of the Cauchy residue theorem has been then extended, and mathematically supported in Ref. [86], where the treatment of dual integrands is modified by the nested application of the residue. Having in mind these features, we constructed a Mathematica package that automates the extraction of these residues by following the lines of Refs. [83,86], the LOop-Tree dualiTY automation (Lotty), and its presentation is the main target of this paper. In effect, several conjectures of results obtained for up-to four loop topologies have always been cross-checked under the umbrella of Lotty.
Furthermore, in a recent work [87], we generalise the findings of Refs. [84,85] and conjectured a close formula for the all-loop causal representation without performing neither an actual application of LTD nor a reconstruction over finite fields. In order to give a support to this conjecture, we implement in Lotty the close formula and give the procedure to generate the causal representation of any loop topology from its features. In fact, we remark that any topology is characterised by the number of cusps and the connections among them, edges. Specifically, the knowledge of causal thresholds can be understood by only considering the structure of the topology, which clearly depends on the number of cusps and edges. This paper is organised as follows. In Sect. 3, we briefly discuss the LTD formalism and set the notation used throughout this work. Then, in Sect. 2, we introduce Lotty and explain how to extract dual integrands, where, to elucidate the discussion, we recall the results for three and four loops found in Refs. [83,85]. We elaborate on the causal representation, in Sect. 4, by presenting the built-in routines in Lotty to symbolically extract causal thresholds and also generate the causal representation of a loop topology from its cusps and edges. As an application, we explicitly show the causal representation of a scalar two-loop double box diagram. Then, in Sect. 5, we discuss the numerical evaluation of causal integrands and present the routines available in Lotty to numerically evaluate them. We remark the integration procedures carried out in Ref. [86] and provide the numerical integration of three-point planar and non-planar scalar triangles at two loops. To check our results, we use the available codes that have implemented the sector decomposition algorithm [99][100][101][102], SecDec 3.0 [50] and Fiesta 4.2 [51]. Finally, in Sect. 6, we draw our conclusions and summarise the features of Lotty.

Loop-tree duality
In this section, we set the notation for the multi-loop Feynman integrals and recap the main features of the novel formulation of the loop-tree duality (LTD) developed in Ref. [83].
In order to start the discussion, let us consider a generic N -point integral at L loops in the Feynman representation, Here and in the following, to simplify the notation in the integration measure, we define Thus, the integral in Eq. (1) is understood as an L-loop topology with r internal lines. Since each internal line is characterised by a propagator of the form q i r = r + k i r (with i r ∈ {1, 2, . . . , r }), where r is the loop momentum identifying this set, and k i r accounts for a linear combination of external momenta { p j } N , one can group lines that share the same integration momenta. For instance, one-, two-and three-loop scattering amplitudes, as shall be described in Sect. 3.4, are characterised by one, three and six combinations of loop momenta, respectively. The numerator N is a function of the loop and external momenta, whose structure is given by the Feynman rules of the theory. The function G F collects all the Feynman propagators of the topology, where α i are the powers of the propagators and, the usual Feynman propagator of a one single particle, with m i its mass, +ı0 the infinitesimal imaginary Feynman prescription and, the on-shell energy of the loop momentum q i written in terms of the spatial components q i . Let us note that in the definition of Feynman propagators (4), we explicitly pulled out the dependence on the energy component of the loop momenta, q i,0 . The latter is carried out in order to profit of the Cauchy residue theorem and integrate out one degree of freedom, which corresponds to the energy component of each loop momentum. Thus, to compute the residues, one selects the poles with negative imaginary part in the complex plane of the energy component of the loop momentum that is integrated, as explained in Ref. [83].
In view of the representation of the Feynman propagators, in terms of the energy component of the loop momenta, the numerator of Eq. (2) can also be expressed as polynomial in these variables, where c are coefficients in terms of the kinematic variables and the spatial components of the loop momenta, α i j are positive integer numbers. This decomposition can be obtained from the d-dimensional scalar product, where q (+) i corresponds to the d-dimensional momentum with on-shell energy q Hence, the dual representation of the integrand (2), after setting on shell the propagators that depend on the loop momentum q i 1 , is defined as follows, where the factor −2πı that comes from the Cauchy residue theorem is absorbed in the definition of the integration measure as shall be noted below. As mentioned before, this residue corresponds to integrating out the energy components of the loop momenta and allows to introduce the nested residue as follows, where the iteration goes until the s-th set and corresponds to setting simultaneously L lines on shell. The latter is equivalent to open the loop topology (or amplitude) into connected trees. Finally, with the integration of the energy component of the loop momentum, one passes from Minkowski to Euclidean space. In the following, we use the abbreviation, for the (d − 1)-momentum integration measure.

Loop-tree duality in Lotty
Lotty -The LOop-Tree dualiTY automation, is a standalone Mathematica package focused on the evaluation of scattering amplitudes by following the ideas of the novel formulation of the loop-tree duality formalism. The routines implemented in Lotty apply, straightforwardly and systematically, the ideas developed in Refs. [83][84][85][86][87]. To illustrate the implementation of LTD described in Sect. 2, we present, in this section, the main functions of Lotty to find the dual representation of multi-loop Feynman integrands together with the decomposition of scattering amplitudes at two and three loops.

Installation
In order to use Lotty, the user needs to have Mathematica installed (version 10 or higher). This package can be downloaded from the public Bitbucket repository, git clone https://wjtorresb@bitbucket.org/wjtorresb/lotty.git Then, to make use of all the features of Lotty, especially for the parallelisation, we recommend defining the variable $LottyPath with the path to Lotty. This can be done by specifying this path in the init.m file of Mathematica. A suggestion for the latter is, Writing the path to Lotty also allows us to simply load this package with the Mathematica command Get["Lotty'"].
In the notebook, one may load the package as follows, Notice that the variables $LottyParallel and $LottyKernels are only needed to parallelise Lotty, otherwise they can be omitted. In fact, by default $LottyParallel = False and $LottyKernels is a variable with an unassigned number. Assuming one is interested in parallelising, the user has to provide the number of kernels on which Lotty is going to be loaded. A naive choice of the latter, using the complete functionality of the computer, is setting $LottyKernels = 2$ProcessorCount.

ResidueW
As briefly recalled in Sect. 2, we need to extensively make use of a routine that analytically extracts the residue of a rational function. Unfortunately, the Mathematica built-in routine Residue turns out not to be efficient when the number of variables of a given rational function starts increasing. Hence, to overcome this issue and make the extraction of residues more efficient, we provide our own version of this routine, ResidueW, which relies on first principles. For instance, the residue of the rational function f (x) with a pole of order m at the point x = x 0 is evaluated as follows, Hence, if f (x) has a single pole at the point x = x 0 , we set m = 1.
Since the denominator of the rational function is constructed from a product of propagators (see Eq. (4)), the extraction of the residue is straightforward due to cancellations that occur before taking the limit of Eq. (11).
Let us illustrate the performance of Residue and ResidueW in a simple rational function with three variables, Clearly, the efficiency of Residue w.r.t. ResidueW starts decreasing when more variables are present in the calculation. For instance, when considering rational functions in terms of more than five variables, e.g. loop topologies with five internal lines, the performance of Residue does not allow to go further in the evaluation of residues. On the contrary, ResidueW overcomes this issue, making the nested application of the residue (9) more efficient and simpler. In all the applications considered so far, dual decompositions of up-to four-loop topologies, we have not observed any limitation in the extraction of the residue.

GetDual
Then, with an improved version of Residue, we can start describing the various routines that automate the LTD formalism. This is the case of GetDual that generates the residue in terms of the on-shell energies, q (+) i,0 , provided that the integrand is a rational function in the energy components of the loop momenta, as displayed in Eq. (6). Hence, the entries of GetDual are the numerator, list of propagators and list of loop momenta. Notice that the explicit dependence of the masses for each propagator is not necessary, since they are included in the definition of q (+) i,0 . To illustrate the use of this function, let us consider a scalar sunrise diagram, without external momenta, at four loops, as the one depicted in Fig. 1  The latter, in the convention of Ref. [83], corresponds to a Maximal loop topology (MLT) configuration at four loops with no insertion of external momenta in the vertices. In general, the MLT configuration is characterised by L + 1 subsets of propagators of the form, where s accounts for the set of all propagators. A further discussion regarding the structure of a loop topology in terms of sets of propagators is provided in Sect. 3.4. Thus, the particular cases of Fig where the dependence on the mass and spatial components of the loop momenta is recovered when expanding the on-shell energies out.
In the former example, we consider, for simplicity, the numerator N = 1, however, the extension to numerators expressed as polynomials in the energy component of the loop momenta, as displayed in Eq. (6), can also be implemented in GetDual. The latter is carried out because the LTD formalism works at the level of energies only. Then, to implement an arbitrary numerator with the explicit dependence on the energy component of the loop momenta, we use the Mathematica function Subscript[#,0]. For instance, the numerators 1,0 and 2 1,0 are simply expressed as Subscript[l1,0] and Subscript[l1,0]ˆ2, respectively. In more details, let us add these numerators to the loop topology described above, Clearly, the former is zero because of the absence of external momenta. The latter, on the contrary, can be understood as, Let us now consider loop topologies with the insertion of external momenta. As mentioned in Sect. 2, the LTD formalism performs the extraction of residues by selecting the poles with negative imaginary part in the complex plane of the energy component of the loop momentum that is integrated, Im q i,0 < 0. Hence, when including external momenta, further assumptions have to be added to account for the latter. Thus, if the loop topology contains external momenta, say p i , one adds for each momentum the assumption Im p i,0 = 0 in the flag "Assumptions" of the routine GetDual. To illustrate this notation, let us consider the same four-loop topology, described above, with the presence of one independent external momentum p 1 , (see Fig. 1 where the dependence on the energy component of the external momentum p 1,0 is completely pulled out after applying LTD. In the applications considered until now, we have dealt with loop integrands in which the dependence on the Feynman propagators is linear. In other words, we have always set α i = 1 in Eq. (2). However, it is pertinent to work with raise powers in the propagators due to the several applications at multi-loop level one finds. To this end, let us consider, yet the same topology, the four-loop integral, where each α i is the i-th power of each propagator. This can be done in GetDual by including the flag "Powers". For instance, the dual representation of the integrand G F 1 2 , 2, 3, 4, 5 can be generated as follows, where, tmp = 1 In Sect. 4, a definition of λ ± i will be provided along the lines of the causal representation of multi-loop Feynman integrands.
Furthermore, we can consider negative powers, which corresponds to dealing with a numerator, e.g.
yielding to,

In[11]:= tmp = GetDual[numerator, propagators, loopmom, "Assumptions" -> assumptions,"Powers" -> {1, 1, 1, 1, -1}]; tmp// Total // ExpandNumerator // FullSimplify
Let us recall that in the applications considered so far, with and without insertion of external momenta, we have summed over all "dual" integrands through the built-in Mathematica routines Total and Together. The explicit structure of individual residues coming from the application of LTD shall be elucidated in the following sections. Besides, these routines in Mathematica become inefficient when increasing the loop order together with the number of internal lines. Hence, as shall be presented in Sect. 4, a representation for the sum of these integrands is given, allowing for straightforward numerical evaluation.

RefineDual
In the previous section, we provide the dual representation of a four-loop sunrise diagram by the direct application of LTD in terms of the on-shell energies q (+) i,0 . However, it would be beneficial for the generation of a dual scattering amplitude to have beforehand a recipe to extract the needed residues of the latter at a given loop order. This study has extensively carried out for topologies up-to three [83,86] and four [85] loops. In this section, we elaborate on these results and show how, within Lotty framework, the parametric form of these residues, regardless of the loop order, can be provided. To this end, we provide the routine RefineDual that performs a "reverse engineering" to find the dual integrands in terms of the various q i,0 .
In order to illustrate the use of RefineDual, let us first recall the LTD formalism at one loop, where the sum runs over all N internal propagators, which are parametrically written as q 1 s = 1 + k s . The latter, in effect, corresponds to considering one set of loop momenta, {q 1 }, with N elements. We define the dual contributions as, where the −2πı factor has been absorbed in the integration measure 1 , previously defined in Eq. (10). The minus sign in −i indicates a reversal of the momentum flow of the corresponding propagators. Notice that this equation is valid independently of the structure of the numerator, provided that the numerator is polynomial in q i,0 . With this in mind, let us recap the most general decomposition of two-and three-loop Feynman integrals in terms of dual contributions, GD.
i. The most general two-loop Feynman integral is characterised by the set of propagators, q 1 s , q 2 s , q 3 s with q 1 s , q 2 s and q 3 s subsets of internal propagators of the form, Then, due to the way how LTD is performed at two loops, two loop momenta are simultaneously set on-shell. Therefore, we can consider, for the sake of simplicity, the subsets q i s with only one element. Thus, by applying LTD through GetDual, 1 that corresponds to the exchange of 1 ↔ 2 w.r.t. the former evaluation. Notably, as shall be discussed in Sect. 4, the ordering in which the residues are computed does not alter the final result, although individual terms are modified. For an extensive discussion on this subject, we refer the reader to Ref. [84].
Although this result is for the particular case of a two-loop topology with three internal lines, it can be generalised to the general case with three subsets, as proven in Ref. [86]. We elucidate this transition in Sect. 3.5 for the two-loop double box topology. A pictorial representation of the lines i is displayed in Fig. 2 (left).
(ii) Similarly, for the most general three-loop Feynman integrals, we have to deal with six subsets, which are pictorially displayed in Fig. 3 (left), whose implementation in Lotty is, The use of this function for any three-loop topology shall be elucidated in Sect. 3.6, by considering the three-loop tennis-court diagram.
The aforementioned topologies, in the convention of Ref [83], correspond to MLT and Next-to-Next-to Maximal loop topology (N 2 MLT) configurations at two and three loops, respectively. Similar to the definition of MLT summarised in Eq. (12), let us briefly recall that the N 2 MLT configuration is characterised by L + 3 subsets of propagators of the form, whose dual decomposition, regardless of the loop order, has been studied by means of convolution and factorisation identities.
In the applications of LTD studied in Refs. [83,85,86], we were interested in achieving a decomposition of topologies of up-to four loops in terms of nested residues. The motivation to do so is to have a parametric representation, regardless of the loop order, in terms of the latter. Thus, allowing to provide the dual representation of any multi-loop scattering amplitude by only applying the corresponding residues. For instance, computing the dual representation of a three-loop scattering amplitude corresponds only to evaluating the residues GD[{i,j,k}] obtained in item ii. We remark that using the routine RefineDual allowed us to cross-check our conjectures and inspired us to provide mathematical proofs for the latter.

Two-loop double box
In the former section, it was discussed that the dual decomposition of any two-loop scattering amplitude can be obtained from the evaluation of a two-loop vacuum diagram (see Fig. 2 (left)), where its internal lines are promoted to be sets of propagators. Hence, to elucidate the connection between the most general loop topology at two loops and the one obtained for a specify two-loop diagram, say two-loop double box, we evaluate the latter as follows, Although we set numerator = 1, the very same analysis can be carried out for any arbitrary polynomial in the energy components of the loop momenta. Thus, evaluating all residues, one gets, Notice that from the list of propagators one can distinguish three subsets according to the dependence on the loop momenta, which corresponds to the same structure of the two-loop residues, parametrised according to Eq. (20) and obtained in item i., after trivially changing α i → i.

Three-loop tennis court
Similar to the two-loop case, let us discuss a particular example at three loops, say three-loop tennis-court diagram (see Fig. 3 (right)), and compare it with the residues obtained from the evaluation of the most general topology at three loops, listed in item ii.
Since the decomposition at three loops is assumed to be valid for any three-loop scattering amplitudes, let us consider, differently from the previous case, the numerator with a i arbitrary constants independent of the on-shell energies. Then, in Lotty, we parametrise the tennis-court diagram as, Based on the organisation of the loop momenta in propagators, we can distinguish the six sets that represent this topology, generating, then, the residues, corresponding to the 12 residues obtained for the three-loop vacuum diagram of Fig. 3 (left) and listed in item ii.

Causal representation in Lotty
In the previous section, we presented the procedure to generate dual integrands in Lotty through the application of LTD.
In the examples considered so far, we have not looked at the structure of each residue in terms of on-shell energies, q (+) i,0 . In fact, we have only been interested in the full sum of these residues. This sum can be easily performed in Mathematica for simple topologies, like MLT configurations (see Fig. 1 for the four-loop case), since it is always expressed as, i,0 ± p 1,0 . The way in which the integrand is written in Eq. (26), in terms of causal propagators, λ ± i , corresponds to the causal representation of multi-loop Feynman integrands. In the recent work of Ref. [87], we provide an interpretation of this pattern by taking into account the features of loop topologies, cusps and edges. To start, we discuss that for the MLT configuration one has two cusps and one edge, where an edge corresponds to a set of all internal lines that connect two cusps, Hence, in the MLT configuration there is only one edge that corresponds to a set with L + 1 elements. Similarly, the N 2 MLT configuration can be seen as a topology made of four cusps and six edges, see Fig. 3 (left). Thus, within this framework, a loop topology made of k + 2 cusps is referred as N k MLT configuration, whose causal representation can be provided along the lines of Ref. [87]. Let us then draw our attention to the N 2 MLT configurations, Fig. 3  Notice that the external momenta satisfy momentum conservation p 1 + p 2 + p 3 + p 4 = 0. Then, the residues are computed as follows, where we introduce the shorthand notation, In order to study the structure of denominators in Eq. (28), let us recall that on-shell energies are positively defined, q (+) i,0 > 0. Therefore, the first and second lines in the residue (28) do not display any unphysical singularity. The only ones that could occur correspond to causal thresholds, given by the values of p i,0 . The third line, on the contrary, displays pseudothresholds or unphysical singularities, which come from the combination of q (+) i,0 's with and without bar, according to Eq. (29), e.g. q (+) (1,4,5),0 = q (+) 1,0 + q (+) 4,0 − q (+) 7,0 . The same pattern is held for all individual residues. Nonetheless, it has been studied in Refs. [84,85] that the sum of all residues completely drops the dependence of all terms that contain at least a bar, q (+) (i 1 ,ī 2 ,...),0 . The former are called non-causal propagators and will not be taken into account in the following discussions.

GetCausalProps
In view of the structure of individual residues, it will be desirable, before summing all contributions up, to know the shape of the denominators in terms of causal propagators, λ ± i = q  [84,85]. In fact, by following the approach explained in [84], together with the analytic reconstruction over finite fields, in the FiniteFlow framework, we manage to reconstruct the causal representation of NMLT and N 2 MLT. Likewise, particular cases of N 3 MLT and N 4 MLT have also been considered in Ref. [85].
In order to provide the causal representation of a given N k MLT configuration, assuming that this topology only has linear Feynman propagators, 2 the integrand in the causal representation admits the following decomposition, Thus, to understand the causal representation, we can only focus on F, since x L+k+1 can straightforwardly be reconstructed at any time of the calculation. In effect, this prefactor can be obtained with the flag "GetPref" of GetCausalProps,

AllCausal
From the former discussion, it is clear that, within LTD, it is cumbersome to have a rational function that displays the structure of integrands and allows for an efficient numerical evaluation. For this reason, our main aim with this causal representation is to find "causal" terms that can straightforwardly be evaluated. On the one hand, there is the subtlety with the unphysical singularities that come from non-causal propagators and, on the other hand, is the evaluation time. To overcome both issues, we only consider expressions written as the one of Eq. (26), where the dependence on causal propagators is explicitly pulled out.
In a recent paper [87], we proposed an algorithm to generate causal representation of topologies constructed from k + 2 cusps. Due to the symmetry in the structure of topologies when all possible connections, edges, between cusps are considered, we conjectured a close formula when all possible edges are considered, 2 Notice that the extension to integrands with raised powers in the propagators can yet be provided by performing, at integrand level, the derivative ∂/∂(q (+) i,0 ) 2 along the lines of Ref. [84]. with, .. contain the causal information of the integrand, The close formula (31)  Notice that while tmp2 contains 24 "causal" terms, tmp3 contains 12. This is because the energies of external momenta are set to zero, corresponding to a vacuum like diagram, whose numerical integration was performed in Ref. [84] in d = 2, 3, 4 space-time dimensions.
In this section, we have discussed the parametric structure of the causal representation of the most general topology in terms of λ ± i without making explicit the dependence on on-shell energies. Clearly, the latter depends on how the propagators for a given topology are defined or, equivalently, how cusps are connected through edges. Let us then elucidate the transition from the all-loop causal representation to a particular one. To this end, we first consider the connection between the four cusps through the six edges (see Fig. 3  where i in q [{i}] corresponds to the i-th cusp in the topology and is defined according to the lines that connect this cusp (red labels in Fig. 3 (left)). Also, in order to pictorially see that the edges defined in SetSingLam are properly written, we implement the routine PlotTop that displays the topology under consideration, producing a figure in a Mathematica notebook, as the one shown in Fig. 4 (left). The assignment of internal lines, according to the cusps, allows to generate the causal propagators λ ± i in terms of q (+) i,0 's. This is internally carried out in Lotty with the function Lamb2qij. For instance, the causal representation of the three-loop N 2 MLT configuration of Fig. 4
On top of the representation found for the three-loop topology, as discussed in Ref. [87], one can profit from the causal representation of the most general topology, in this case N 2 MLT, and obtaining without any evaluation of the residue the causal representation of topologies with lower number of edges, but same number of cusps. In effect, we can also generate topologies with five ( Fig. 4 (center)) and four (Fig. 4 (right) Then, to obtain the causal representation of a two-loop kite diagram (Fig. 4 (middle)) with two external momenta, one focuses only on the connection of edges and momentum conservation, Let us remark that, as explained in Ref. [87], cancellations due to the absence of q (+) 5,0 and/or q (+) 6,0 appear. In particular, in the former cases the number of "causal terms" decreases from 24 to 20.

RefineCausal
It is clear that when generating the most general causal representation of a loop topology, given by a fixed number of cusps, several redundant terms will appear. Hence, to compensate their appearance, we work out the expressions of L ± before expanding in λ ± 's. This is straightforwardly carried out by studying the connections between cusps, provided e.g. by SetSingLam2 and SetSingLam3, through the routine RefineCausal.
To illustrate the use of RefineCausal, let us keep studying the topologies with four cusps of Fig. 4, where one-and two-loop topologies are particular cases of the N 2 MLT configuration. Then, starting from the general causal representation, we explicitly add the connections between cusps and edges for the two-loop topology as follows, This result is expected since the cusps 1 and 4 do not share any edge, as noticed in the definition of SetSingLam2, and allows for reducing from 24 to 22 causal terms. Besides, in this configuration, due mainly to momentum conservation of the external momenta, we have a further relation, by taking into account the explicit dependence on λ ± in Eq. (33) and λ ± 14 = λ ∓ 23 . Thus, recovering the 20 causal terms of Ref. [87].
Similar for the one-loop box, representation proposed in Ref. [87]. Clearly, the former will display several numerical instabilities w.r.t. the latter. Then to avoid this issue, we evaluate both integrands in the non-physical region, Im( p i ) = 0. In order to compare both approaches, let us focus on the one-loop box of Fig. 4 (right), and generate the dual representation of this integrand, Notice that, in order to find a numerical agreement between causal and non-causal (sum of residues) representation, we profit from the evaluation of the integrand over random rational points. The same approach is carried out for all the examples collected in usage.nb.

Six-cusp topologies
In the previous section, we discussed that the routine AllCausal allows to generate the most general causal representation of a loop topology, given by the number of cusps and their connections through edges. Remarkably, the use of LTD is not needed to generate those expressions, but it is important to numerically compare both approaches. Hence, to elucidate this advantage, in this section we present the causal representation of loop topologies made of six cusps. The general causal representation of a topology with six cusps generates topologies up-to 15 loops. In Fig. 5, we display a few topologies (up-to three loops) that can be obtained from the most general one. In the following, we elaborate on the two-loop double box. The other topologies at one and three loops are analysed in the notebook usage.nb, following the same approach as this section.
Let us then start with the list of propagators of the two-loop double box of Fig. 2 (right) In [42] := propagators = {l1, l1+p[1

], l1+p[1]+p[2], l2, l2+p[4], l2+p[3]+p[4], -l1-l2};
where, by means of cusps and edges, this topology can be understood as, The edges in SetSingLam are defined according to Fig. 2 (right)  Similar to the one-and two-loop topologies considered in the previous section, in tmp3, we generate the most general causal representation for a topology built from six cusps. Then, in momcon, we define the momentum conservation of the four external particles, setting to zero those that are not present in this configuration, p 5 and p 6 . The routine Lamb2qij, as explained above, performs the replacement from λ i 's to q (+) (i 1 ,i 2 ,...),0 . Finally, in the last line, we apply all the conditions. Since the full value of tmp4 is not illuminating for the current discussion, we just verify the variable the latter depends on, In order to check that the stored expression in tmp4 corresponds to a better rewriting of the dual integrands in terms of causal propagators, we compare this result with the one obtained directly from the nested application of the residue (9). To this end, we generate the residues in the standard way, as done in Sect. 3.5, whose difference in both evaluations is an overall sign, which can be accounted for in the general expression of Eq. (31).

Numerical integration
In the previous sections, we discussed and provided automated tools to generate integrands in the causal representation. As mentioned throughout this paper and in former works, the aim of doing so is to be able to perform smooth and efficient numerical evaluations of the latter. Then, with this approach in mind, we provide non-trivial examples of the evaluation of three-point finite integrals at two loops. Differently from previous studies, we consider multi-loop Feynman integrands with presence of external momenta, yet getting agreement with publicly available softwares based on sector decomposition. Since the main feature of Lotty is performing as many symbolic operations as possible to render in a very simple integrand, we implement, for numerical purposes, the change of variables discussed in Ref. [84], which we recall presently.

Change of variables
In order to perform a numerical integration of an L-loop Feynman integrand, one has to deal, in general, with (d − 1)L integrations, which can be straightforwardly embedded in (d − 1)L-dimensional hypersphere or in the product of L (d − 1)dimensional spheres. 3 The former and latter change of variables are implemented in Lotty by activating the flag "Domain" with the options "HyperSphere" and "Sphere", respectively, in the routines LoopToSC, meassure and intvars.  Let us comment on the routines we introduced in the previous example. LoopToSC provides the replacement rules to spherical coordinates for the spatial components of the loop momenta as (d − 1) arrays for each loop. meassure and intvars write, according to the chosen "Domain", the integration measure and integration variables, respectively. Notice that, due to the LTD formalism, the simplicity in this change of variables stems of the transition from Minkowski to Euclidean space. The default domain of the above routines is "Sphere", which, as mentioned before, corresponds to embed all integrations in a product of L (d − 1)-dimensional spheres. In this case, the notation for angles θ i j and radius r i exactly corresponds to the i-th loop momentum in the list LoopMom. In the integration limits of each radius, one can choose the upper bound, by default we set "MaxValueR"− > 10 5 . Nevertheless, to avoid considering the integration domain [0, ∞), one can map [0, ∞) to the segment (0, 1], with the change of variables, making thus the numerical evaluation of the integrand simpler. Since the integrands are written in terms of the on-shell energies q (+) i,0 , we provide an additional routine to perform the scalar products in the Euclidean space, sp [a,b]. For instance, the scalar product of 1 · 2 in (4 − 1) dimensions, where, to evaluate these scalar products, one needs to add two lists (with the same length). In the former example, threedimensional arrays for 1 and 2 are given by the parametrisation, in spherical coordinates, provided by LoopToSC[l1, l2, dim], as mentioned above. Hence, with the routine sp ready to perform (d − 1)-dimensional scalar products in Euclidean space, we provide, for illustrative reasons, the on-shell energies of the two-loop planar triangle of Fig. 6  Notice that the user has to include the four-dimensional arrays for the external momenta. In this particular example, the external momenta are chosen to be in the center-of-mass frame (COM). Thus, Ecm corresponds to the COM energy, yielding to the on-shell energies in spherical coordinates,  where in the temporary variable value we make the replacement assuming that the mass of the i-th propagator is m i . In the derivation of the on-shell energies, let us note that one can harness from the choice of the reference frame. For instance, in the evaluation of the two-loop planar triangle of Fig. 6 (left), one has that the only massless propagator is 1 + 2 while the other ones are massive with the same mass. Then, working in the COM frame, a simplification at the level of on-shell energies is manifest, q 4,0 . This choice of the reference frame reduces the number of variables to evaluate from 6 to 4. In the following section, we discuss the stability of the causal representation of this integrand.

Numerical stability
Then, with q (+) i,0 expressed in spherical coordinates, we can study the numerical stability of integrands in the causal representation. To do so, we consider the two-loop planar and non-planar scalar triangles depicted in Fig. 6, whose complete study of the renormalised amplitudes have been carried out in Refs. [68,69] for the Higgs production in QED and in QCD/EW, respectively.
In the following examples, we work in the region where no threshold singularities appear in the evaluation of the integrand. Nevertheless, the extension to integrands containing thresholds can yet be carried out with the aid of the routine GetCausalProps presented in Sect. 4, since the latter explicitly displays the causal propagators and thus the locus of all thresholds. Then, we specify the mass of propagators according to Fig. 6, where, for the planar triangle, the propagator with a curly line is considered as massless whereas the other ones are massive. For the non-planar triangle, we choose all masses to be the same. We recall that choosing different values of masses for the internal lines (solid and wavy) do not generate any limitation to our approach, since we perform a purely numerical evaluation, as observed in Ref. [69].
In Sect. 5.1, we provide the change of variables to spherical coordinates leaving the values of radii r i still unconstrained and also discuss the mapping to x i , [0, ∞) → (0, 1], to avoid the boundary limits. With this in mind, we produce, in Fig. 7, a 3D plot to show the numerical behaviour of the planar and non-planar two-loop triangle integrands. We vary simultaneously x 1 and x 2 , keeping the angles θ i j fixed. Differently to what is carried out in former works [84,85], we are not only focused on the structure of the integrand in terms of q (+) i,0 . Instead, we provide the complete behaviour of the integrands, where because of the absence of UV and IR singularities in the latter, a numerical integration can straightforwardly performed. On top of the 3D plots, we also provide plots in 2D, in Fig. 8, where we scan over x 1 (x 2 ) keeping fixed x 2 (x 1 ) and varying the angle θ 12 , yet displaying a smooth behaviour at integrand level. In other words, the 2D plots in Fig. 8 correspond to slices of Fig. 7 for a given x i .
At this level, we are left with the numerical integration. Then, with the Mathematica built-in function NIntegrate, we present a coarse integration in which neither an optimisation nor an educated use of the latter function was carried out. The results of numerical integrations of several phase-space points, compared with SecDec 3.0, are collected in Table 1, where agreement up-to three digits is found with a simple numerical evaluation. The evaluation time per point is O 30 in a laptop machine with an Intel i7 (1.2 GHz) processor with 4 cores and 8 GB of RAM. We provide a Mathematica notebook, LTD_integrations.wl, that can be downloaded together with Lotty, with the complete calculation presented in this section.
We remark that the main purpose of Lotty is not providing a Monte Carlo to numerically integrate our integrands. On the contrary, the aim of this Mathematica package is to give a straightforward understanding of the LTD formalism and,  Fig. 6 with random samplings of x 2 and x 1 in the left and right plots, respectively therefore, the causal representation of multi-loop Feynman integrands. In view of the latter, and the various routines to compute the dual integrands by means of the Cauchy residue theorem, we claim that any numerator written as polynomial in the energy component of the loop momenta can be treated along the lines discussed in this paper. Thus, promoting the causal representation of integrands to scattering amplitudes. Moreover, it will be desirable to understand alternative strategies that allow for a more stable numerical integration.

Summary and outlook
The novel formulation of the loop-tree duality (LTD) allows for an alternative approach to numerically evaluate multi-loop Feynman integrals. Interestingly, the representation of dual integrands by means of LTD has been found to display only physical information at integrand level. Therefore, the appearance of non-causal contributions, or often called pseudo-thresholds, does not generate obstacles when numerically evaluating the integrands generated through LTD.
In view of this novel formulation and the results that are being obtained in Refs. [83][84][85][86], for topologies up-to four loops and their generalisation to L loops, it is desirable to provide a stand-alone package that performs all the calculations summarised in former works. To this end, in this paper we presented the Mathematica package Lotty that automates the novel formulation of LTD and elaborates on the causal representation of loop topologies. In particular, the close formulae for the all-loop causal representation conjectured in Ref. [87] were implemented in Lotty, giving, in this way, a support to the results presented in the latter.
In order to make a connection between the formalism proposed in Ref. [83] and the automation provided by Lotty, we concentrated on particular cases at two and three loops, showing that the conjectures initially proposed in Ref. [83] and then proven in Ref. [86] can be worked out regardless of the loop order. In particular, we discussed the dual representation of two-and three-loop scattering amplitudes, showing explicitly that this decomposition, obtained for vacuum-like diagrams, is generalised to topologies with higher multiplicity in the external momenta. In order to illustrate this behaviour, we considered the two-loop double box and the three-loop tennis-court diagrams.
Elaborating on the all-loop causal representation conjectured in Ref. [87], we presented the routines in Lotty that provide a parametric form of the integrands in terms of the physical singularities of the topology. To ensure that the latter formula is correct, we numerically compared the results obtained from the direct application of LTD against the causal formulae, finding completely agreement. This numerical evaluation was performed by considering samplings of random rational numbers, in order not to lose any precision in the evaluation.   Furthermore, to study the numerical stability of the integrands in the causal representation, we considered, as examples, the evaluation of two-loop planar and non-planar triangles, in which the change of variables to spherical coordinates, implemented in Lotty, allowed us to study the numerical stability of these integrands. In particular, we presented two-and three-dimensional plots displaying their numerical stability.
Then, in view of the smooth behaviour of these integrands, we numerically integrated out our expressions in several phase-space points. This integration was naively performed in Mathematica through the function NIntegrate, finding agreement with the results provided by the automated code of SecDec 3.0.
The automation presented in this paper allows to implement in a simple way the LTD formulation and generate multi-loop integrands containing only physical information. In effect, the dual decomposition of topologies up-to four loops, including presence of external momenta, have been tested without showing a limitation. Likewise, in the generation of the all-loop causal representation of a topology given by its features, cusps and edges, Lotty has been able to provide causal representation of loop topologies made of up-to fifteen cusps.
In view of the implementation of LTD and the structure of the causal integrands generated by Lotty, one can focus only on dealing with physical singularities, infrared (IR), ultraviolet (UV) and thresholds, and aim at a four-dimensional finite object at integrand level. The cancellation of these singularities has been successfully performed at next-to-leading order (NLO) through the four-dimensional unsubtraction scheme [70][71][72] and it is under investigation at next-to-next-to-leading order (NNLO). In particular, the calculation of IR-safe objects at NNLO started to appear in the numerical evaluation of three-point scattering amplitudes at two loops [68,69]. More applications, including the treatment of IR singularities, have to be considered to provide a purely numerical implementation at NNLO.
We expect that the use of LTD together with its automation, given in this paper by Lotty, pushes forward the calculations of physical observables that currently are displayed obstacles when trying to compute them with the standard methods. We emphasis that Lotty is ready to be interfaced with publicly available softwares that perform the generation and calculation of multi-loop Feynman integrands, e.g. Refs. [103][104][105][106][107][108].
Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: Our manuscript does not have external data associated to it].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .