Numerical Loop-Tree Duality: contour deformation and subtraction

We introduce a novel construction of a contour deformation within the framework of Loop-Tree Duality for the numerical computation of loop integrals featuring threshold singularities in momentum space. The functional form of our contour deformation automatically satisfies all constraints without the need for fine-tuning. We demonstrate that our construction is systematic and efficient by applying it to more than 100 examples of finite scalar integrals featuring up to six loops. We also showcase a first step towards handling non-integrable singularities by applying our work to one-loop infrared divergent scalar integrals and to the one-loop amplitude for the ordered production of two and three photons. This requires the combination of our contour deformation with local counterterms that regulate soft, collinear and ultraviolet divergences. This work is an important step towards computing higher-order corrections to relevant scattering cross-sections in a fully numerical fashion.


Introduction
The Large Hadron Collider (LHC) is entering its high luminosity data acquisition phase and is thus transitioning from being a discovery experiment to a precision measurement one. For this new goal, accurate theoretical predictions are necessary in order to ensure that theoretical uncertainties remain at or below the level of experimental ones. In particular, this involves the computation of higher-order corrections to the cross-sections of relevant scattering processes, which are built by considering processes with additional unresolved partons (real-emission type of contributions) and additional loop degrees of freedom (virtual type of contributions). These two classes of contributions are separately divergent but combine into a finite quantity in virtue of the Kinoshita-Lee-Nauenberg theorem [1,2]. Traditionally, the computation of these two components is performed using very different approaches and the deep connection relating their degenerate infrared degrees of freedom is only realised through dimensional regularisation [3][4][5] at the very end of the computation. Indeed, real-emission contributions are typically computed numerically through the introduction of subtraction counterterms [6][7][8][9][10][11][12][13][14][15][16][17] or some form of phase-space slicing [18][19][20][21][22][23][24][25], whereas the evaluation of their virtual counterparts is mostly carried out purely analytically, thus realising the cancellation of infrared singularities at the integrated level. A notable exception is the computation of inclusive Higgs production at N 3 LO accuracy [26], which was performed through reverse-unitarity [27,28]. Even though the same technique was further developed to accommodate the Higgs rapidity distribution in ref. [29], it is clear that this approach is not applicable to fully differential high-multiplicity processes. Furthermore, despite impressive advances in the mathematical aspects of the reduction of scattering amplitudes to master integrals [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45], and their subsequent computation by means of differential equations [32,[46][47][48][49][50][51], it is believed that the computation of many relevant higher-order corrections to important processes (e.g. NNLO corrections to pp → ttH and pp → ttbb) will remain intractable with this traditional approach, in part due to the increase in the number of scales relevant to the problem and because of the appearance of new mathematical structures in the form of generalised elliptic polylogarithms [52][53][54][55][56].
Numerical alternatives have been developed for the direct evaluation of loop integrals through sector decomposition [57][58][59][60][61][62] of their Feynman parametrisation or semi-numerical solutions [63][64][65] of the system of differential equations relating them. This lead to the flagship computations of the NNLO corrections to the processes pp → HH [66,67] and pp → tt [68], where the exact dependency on all quark masses was kept. Although these achievements demonstrate the superiority of numerical approaches in selected cases, they still suffer from the scalability issue inherited from their reliance on the analytical reduction of the complete amplitude to master integrals.
In light of the above overview of the research field of precise collider predictions, we choose to pursue an alternative construction which considers a purely numerical integration of the virtual contribution in momentum space. One particular benefit from such an approach lies in the prospect of bypassing the reduction to scalar integrals by considering the numerical integration of complete amplitudes directly (see existing results for one-loop amplitudes in refs. [69][70][71][72] and first steps for applications to higher-loop finite scalar integrals in ref. [73]). Working in momentum space is especially appealing when also performing the loop energy integral(s) analytically using residue theorem. This energy integration yields the Loop-Tree Duality (LTD) which provides an alternative representation for the loop integral containing terms with as many on-shell constraints as there are loops, making them effectively trees. This aligns the measure of phase-space and LTD integrals, thus making LTD ideally suited to pursue the ambitious goal of directly combining real-emission and virtual contributions and compute them numerically at once by realising the local cancellation of their infrared singularities. As with reverse-unitarity, this direct-unitarity treatment explicitly maintains the aforementioned connection between real-emission and virtual contributions which is lost when computing them separately or using Feynman parametrisation. Pioneering work of ref. [74][75][76][77][78] demonstrated the potential of carrying out this numerical programme by applying it at one loop. However, during the last decade, the NLO revolution and the successes of analytical methods for the computation of many NNLO-accurate 2 → 2 cross-sections mostly overshadowed such purely numerical approaches. That is until recently, when groundbreaking new results from traditional analytical techniques arguably slowed down, thus opening the way for more numerical alternatives.
Since such radically different purely numerical approaches have to be developed from the ground up, they will not immediately catch up with the impressive analytical work performed by the community over the last two decades. Instead, we proceed incrementally and build progressively towards the complete numerical evaluation of higher-order corrections while making sure at every step that our partial results are robust and make no compromise regarding generality in terms of the perturbative order and process considered. We started this endeavour with ref. [79] where we derived a general formulation of LTD by iteratively applying one-dimensional residue theorem. We showed how the duality relation hence obtained can easily be constructed algorithmically for any loop count and topology, and we tested it by applying it to many integrals without threshold singularities. In that regime, we could perform the integration of the LTD integrand directly as it does not require any contour deformation or counterterms.
The first part of this work concerns the natural follow-up to ref. [79]: regulating threshold singularities in order numerically integrate loop integrals evaluated with physical kinematics. We achieve this by constructing a contour deformation in the (3n)-dimensional complex integration space, designed in accordance with the constraints imposed by the causal prescription of Feynman propagators and by the matching conditions stemming from analytic continuation. Contour deformations for numerical integration have been considered in the past [69,73,80], and we present a novel variant well-suited to our multi-loop LTD expression. In order to ensure that our construction is correct for arbitrary (multi-)loop integrals, we apply it to more than a hundred qualitatively different examples, always finding agreement with the analytical benchmark (when available). We also demonstrate in this way that the convergence rate of our current numerical implementation already renders it competitive. Finally, we discuss optimisation strategies to explore in future work that can improve results further.
The second part of the paper is dedicated towards applying our numerical programme to the computation of divergent scalar diagrams and of physical amplitudes. We consider divergent scalar box and pentagon topologies and the one-loop correction to the ordered production of two and three photons from a quark line. This amplitude involves soft and collinear singularities that correspond to pinched threshold singularities where no regulating contour deformation is allowed. This type of singularities can therefore only be regulated by the introduction of ad-hoc counterterms or through a direct combination with real-emission contributions. In this work, we consider the former. In the case of scalar integrals we introduce a method to remove all IR divergences from one-loop diagrams. When considering complete amplitudes we combine our contour deformation and LTD integrand with the infrared and ultraviolet counterterms presented in ref. [81,82].
The outline of this work is as follows. In sect. 2, we fix our notation by recalling our general multi-loop LTD expression. We construct a general contour deformation in sect. 3. In sect. 4, the subtraction procedure for one-loop scalar integrals and amplitudes is discussed. In sect. 5, we discuss various optimisations for our numerical integration. In sect. 6, we discuss our numerical implementation and we show our results in sect. 7. Finally, we present our conclusion in sect. 8.

Loop-Tree Duality
In this section, we fix the notation and summarise our findings presented in ref. [79]. A general n-loop integral in four-momentum Minkowskian space can be rewritten as an integral over the Euclidean space of the three-dimensional spatial part of the loop momenta. The integrand in that case is the sum of residues obtained by iteratively integrating out the energy variables one after the other by applying residue theorem. Each residue identified in this manner corresponds to a particular spanning tree (i.e. a tree graph that connects all vertices) of the underlying loop graph, or equivalently, to a particular loop momentum basis (i.e. the n edges that complete a spanning tree back to the original n-loop graph) together with a specific set of signs for the energy solutions of the on-shell conditions fixing the residue location, which we call the cut structure.
More precisely, we start from the following n-loop integral where e is the set of indices labelling the edges of the connected graph identifying the integral considered and the numerator N is a regular function of the loop momenta. We assume the Feynman propagators to be pairwise distinct with on-shell energies ±E i = ± q 2 i + m 2 i − iδ. The momentum flow in a graph is uniquely determined by the choice of (consistent) signature vectors s i = (s i1 , . . . , s in ), s ij ∈ {±1, 0} for each propagator, such that q µ i = n j=1 s ij k µ j + p µ i , where p µ i is a shift that depends on external momenta. We consider the integration of the energies in a fixed arbitrary order, set by (k 0 1 , . . . , k 0 n ), each along the real line 1 and closing on an arc of infinite radius in either the upper (with winding number Γ j = +1) or the lower (Γ j = −1) complex half-plane. We assume the integrand to vanish for large loop momenta, so that we can consider the integral along this arc to be zero, thus allowing us to relate the original integral to the sum of residues at poles located within the contour.
When carrying out this iterative integration of the loop energies and collecting residues, one finds that some residues may lie within or outside the integration contour depending on the spatial part of the loop momenta. This would be an unfortunate complication, but we conjectured and verified explicitly that only the residues that unconditionally lie within the integration contour contribute to the integral, and moreover with the same prefactor, whereas all other conditional residues are subject to exact cancellations [79]. We write the dual integrand corresponding to one particular residue of the original integrand f = N/ i∈e D i identified by the loop momentum basis choice b = (b 1 , . . . , b n ), b j ∈ e (corresponding to the list of propagators put on-shell for this residue) as It describes a residue that is within the contour for all loop momentum configurations if where

4)
1 As discussed in ref. [83], our final expression in eq. (2.5) is also correct in the case of complex-valued external momenta, due to the fact that the right-most column of the matrix appearing in eq. (2.4) does not include the imaginary part Im[p 0 i ] of the external momenta. We note however, that the correct interpretation of the absence of this term in eq. (2.4) for complex-valued external kinematics is that the energy integrals are no longer performed along the real line but instead along a path including only one out of the two complex energy solutions of each propagator.
which for a choice of integration order, contour closure and momentum routing (determined by ( k 1 , . . . , k n ), Γ j and s ij respectively), is satisfied unconditionally for exactly one configuration of signs, the cut structure, denoted by σ b . Therefore, the original integral of eq. (2.1) is identically equal to the resulting LTD expression (2.5) where B is the set of all loop momentum bases. We stress again that the functional form of the LTD expression is implicitly dependent on the chosen order for the integration of loop energies, the contour closure choices and the particular momentum routing chosen for the original integral. However, we verified explicitly that one always numerically obtains the same result for the sum of residues for given values of the spatial part of the loop momenta (set in a particular basis). In order to facilitate the understanding of the central result of eq. (2.5), as well as to give some insight on its derivation, we provide an explicit two-loop example in appendix A. Finally, we provided as ancillary material of ref. [79] a Python implementation of the automated derivation of the cut structure for arbitrary loop topologies. Beyond its practical value, this code also demonstrates that explicitly unfolding eq. (2.5) can be done without any computational overhead.
The dual integrands can become singular on surfaces which may be labelled by the residue corresponding to the particular dual integrand in which they appear (specified through the loop basis b) and the particular propagator of that dual integrand that becomes on-shell (specified through the propagator index i). These singular surfaces are of the form i are implicitly defined through the change of basis q µ i = j∈b s b ij q µ j +p µ,b i induced by the loop momentum basis b identifying this surface. The singular surfaces ξ can be separated into two classes: E -and H-surfaces. E-surfaces are defined by the property of having all signs α k , k ∈ b ∪ {i} equal, unless α k is zero. We call the particular sign that all α k are equal to (when not being zero) the surface sign. We factor out the surface sign and name the resulting E-surface η b,i . From this point on, we consider every E-surface to have a positive sign for all energies: E-surfaces are convex and bounded. H-surfaces are then defined by having at least one positive and at least one negative α k and they are labelled γ b,i,α i . A particularly elegant feature of LTD is that the sum of dual integrands forming eq. (2.5) only becomes singular on E-surfaces, as the singularities from H-surfaces cancel pairwise thanks to a mechanism referred to as dual cancellations [84,85]. For δ = 0, an -5 -E-surface has a non-empty set of real solutions in k = ( k 1 , . . . , k n ) ∈ R 3n if it satisfies When both sides of this inequality are exactly zero, the E-surface has no interior since its minor axis is zero, and the E-surface corresponds to the location on an infrared collinear and/or soft singularities of the integral. We refer to them as pinched E-surface, with the important property that singularities they correspond to cannot be regularised via a contour deformation of the loop momenta integration phase-space. For δ > 0 an E-surface η is uniquely regulated by the imaginary prescription We do not find it particularly useful to work out the imaginary part of the the squared propagators appearing in eq. (2.5) (referred to as dual propagator in ref. [84]). Instead, we prefer to stress that the relevant imaginary part of the E-surface equations induced by the causal prescription has a simple definite sign. As it will be made clear later, this observation is indeed the only relevant one in regard to the construction of a contour deformation that satisfies physical requirements and regulates threshold singularities.

Contour Deformation
Numerical integration of Feynman diagrams and physical amplitudes in momentum space originated with the early attempts by Davison E. Soper in [86] and [87], in which the LTD formalism was applied to virtual diagrams at one loop in order to then integrate the cross-section directly. Interestingly, the author also explicitly mentions and utilises the mechanism of local real-virtual cancellations to render the integrand finite at the location of the non-integrable soft and collinear singularities. In order to avoid so-called scattering singularities, referred to in our work as one-loop E-surfaces, the author devised a contour deformation capable of satisfying the relevant constraints. Several methods have since been developed for integrating diagrams and amplitudes directly in four-dimensional loop momentum space. A first success was the computation of one-loop photon amplitudes in ref. [88], followed by refs. [70,[89][90][91] which generalised the formalism beyond one loop and applied it to more challenging integrals. The especially inspiring feature of this series of publication is the focus on constructing a provably exact deformation, through the concept of anti-selection and dynamic scaling of the deformation.
Around the same time when these techniques were developed, a different line of work expanded on LTD and, specifically, on its aspects relevant for the (3n)-dimensional numerical integration of integrals, amplitudes and cross sections [92][93][94]. The contour deformation presented in these works is based on a linear combination of vectors normal to the existing E-surfaces, weighted by adjustable parameters and dampened by exponential functions with unspecified width; the deformation proves to be correct for simple threshold structures and in the limit of arbitrarily small dampening widths. Results obtained in this way however -6 -highlighted for the first time the potential of numerical integration over the spatial degrees of freedom resulting from the LTD identity.
In this section we will construct a reliable and exact deformation that is valid for an arbitrary number of loops and legs. We will give specific examples in order to illustrate how to implement the deformation constraints for complicated singular structures, especially on intersections of multiple E-surfaces.
As long as an integral only features non-pinched threshold singularities, it is possible to engineer a contour deformation yielding a finite result for the integral. The absorptive part of the integral is correct provided that the contour deformation considered satisfies requirements imposed by physical conditions, in particular causality. In relativistic quantum mechanics, causality is originally realised in Feynman propagators via the iδ-prescription or, equivalently, by the request that the theory is in the range of validity of Gell-Mann and Low's theorem [95]. In the LTD formalism, an imaginary prescription on propagators remains and, although its formal expression is more complicated than iδ, it still holds that on E-surfaces this prescription sign is fixed (i.e. it does not depend on either external nor loop kinematics, see eq. (2.9)).
Contour integration of threshold singularities requires to analytically continue the LTD integrand by replacing its dependence on the chosen basis of loop momenta k, by the com- The spatial momenta associated with each propagator are a linear combination of the vectors in the chosen loop momentum basis plus an affine term: Once analytically continued, these spatial momenta then also acquire an imaginary part: Each surface η has an associated energy shift p 0 η , defined in eq. (2.7) as a specific linear combination of the energies of external particles.
An approximation of the imaginary part of the E-surface η can be obtained from the first order term of its Taylor expansion in κ : on E-surfaces (eq. (2.9)) with the sign of the imaginary part of E-surfaces that results from the deformation (eq. (3.3)). We obtain: Magnitude: The norm of the deformation vector is limited by three constraints: Integrand continuity: The LTD expression can be seen as a function of the on-shell energies of the internal particles E i = q 2 i + m 2 i − iδ. These square roots have to be evaluated on a well-defined Riemann sheet. Thus the contour must not cross the branch cuts of any of the involved square roots.
Complex pole constraint: By extending the domain of the LTD integrand from R 3n to C 3n through the replacement of its functional dependency on k with ( k, κ), we find that in addition to real-valued poles (corresponding to the existing E-surfaces), the integrand also features complex-valued poles located at ( k, κ), with κ = 0. We stress that these complex poles exist for all E-surface equations: those (pinched or not) already having solutions for real loop momenta ( k, 0) as well as those that do not and which are referred to as non-existing E-surfaces (in regard to the fact that their existence condition of eq. (2.8) is not fulfilled).
According to Cauchy's theorem, the result of the contour-deformed integral will only be identical to that of the original defining integral over the spatial part of the loop momenta in the real hyper-plane, if and only if the volume defined by this real hyper-plane and the deformed contour does not contain any of such complex poles. The magnitude of the contour deformation must therefore be constrained to be small enough so as to exclude these complex poles.
Expansion validity: The causal constraint on the direction of the contour deformation as well as the complex pole constraint are derived from the Taylor expansion of each energy function E i . We must therefore impose that the norm of the contour deformation vector field is such that the complex argument of each square root defining an energy remains within the range of validity of its expansion.
The next section 3.1 presents the one-loop contour deformation direction constraints and our approach for solving them. We will refer explicitly to illustrative examples that introduce key concepts of our work. The precise and complete description of our construction of a contour deformation valid for an arbitrary number of loops and legs is presented in sect. 3.2.

Pedagogical construction at one loop
Consider a one-loop scalar box diagram in the LTD representation after having explicitly solved the on-shell constraint: where we used that at one loop the dual propagator factorises into the product of an Eand an H-surface, as At one loop, one can also simplify the loop basis identifier b and write it as the index b ∈ e = {1, 2, 3, 4} corresponding to the single LTD cut considered. Thanks to the mechanism of dual cancellations, the sum of all dual integrands is only singular on E-surfaces which, at one loop, are two-dimensional rotational ellipsoids in spatial loop momentum space. All of the potential singular E-surfaces of this scalar box appear as zeros of the functions with i, b ∈ e, i = b, and for given four-momenta of the four external legs p ext j , j ∈ {1, 2, 3, 4}. The number of E-surfaces that have solutions for real loop momenta has an upper bound based on the topology and the number of legs N . For one-loop topologies, an upper bound on the total number of existing E-surfaces is N (N − 1)/2, since we require b = i and using the fact that if η bi exists, η ib cannot exist.
The singularity structure of the LTD expression can be studied by focusing on particular singular E-surfaces and their intersections. In order to do this, we define the boundary and interior operators as The E-surface η bi exists, that is If two ellipsoids η, η exist and intersect, then ∂η ∩ ∂η = ∅. Furthermore, if they intersect without being tangent, they also overlap: ∂ − η ∩ ∂ − η = ∅. As an illustrative example, we now set particular values for the external box kinematics, which we refer to as and list the resulting four members of the set of existing E-surfaces E = {η 12 , η 13 , η 42 , η 43 }, It has four singular E-surfaces with four (partially shared) focal points coloured in red.
The four E-surfaces in eq. (3.11) are coloured according to the colour scheme used in fig. 1. A focal point is the loop momentum (k x , k y , k z ) that sets the argument of an energy square root to zero. Each ellipsoid has two focal points, indicated with red dots in the figure.
The energy shift p 0 i − p 0 b is the length of the major axis. The particular external kinematic configuration chosen in eq. (3.10) has no component along the k z -axis and therefore the particular section k z = 0 corresponds to the plane where the four E-surfaces have a maximal extent.
According to eq. (2.9) we require the imaginary part on any E-surface η to always be negative: sgn(Im[η]) = −1. By replacing k → k − i κ( k) and expanding the E-surface equations to first order in || κ||, we find that the prescription reads which imposes that on any point on the E-surface, κ( k) should point outwards of the Esurface. On the intersection of many E-surfaces, the combined prescriptions impose that κ( k) must simultaneously point outwards of all of the intersecting E-surfaces. One choice that always satisfies the condition of eq. (3.12) for one single E-surface as well as for two intersecting E-surfaces is the sum of their respective normal vector fields, as shown in fig. 2. A similar deformation was proposed in ref. [93], where the deformation field κ( k) is written as a linear combination of the normal fields weighted by an exponential dampening factor that ensures that each normal field vanishes away from its defining Esurface. This particular choice of deformation vector is unsatisfactory when more than two -10 - In fig. 3 we give an example with three E-surfaces, where a naive unweighted sum of normal vectors does not yield a valid deformation. By using fine-tuned dampening of the normal vector fields from each E-surface, such cases may be avoided but this does require an ad-hoc treatment and can lead to poor numerical convergence.
The next subsection introduces the concept of deformation sources which we will use to build a deformation that avoids the shortcomings discussed in this section when considering normal fields.

Deformation sources
Since E-surfaces are convex surfaces, given a point s within the interior of an E-surface ∂ − η, the radial field v s ( k) ≡ k − s, centered at s, satisfies the causal prescription Im[η]| k−i v s < 0 on any point on the surface, where η( k) = 0. We note that the interior of the intersection of a set F ⊆ E of E-surfaces again defines a convex volume and therefore we analogously have that, for any given point s in this volume, that is s ∈ η∈F ∂ − η, the corresponding radial field v s simultaneously satisfies the causal prescription of all on the E-surfaces in F and, especially, on their intersections. We call such a point s a deformation source for the overlapping set F . For a case in which there exists a single point s simultaneously in the interior of all of the existing E-surfaces, then the radial deformation field κ( k) ∝ ( k − s) satisfies the causal prescription on all the threshold singularities.
When there is no single point simultaneously in the interior of all E-surfaces, one can construct a deformation vector written as the sum of radial fields centered at different locations, and adequately multiplied by an anti-selector function disabling the effect of the radial field on all the E-surfaces in which the point is not contained. The anti-selection is constructed such that the individual terms building the deformation vector fields are always "additive" in their ability to satisify the causality requirements. Indeed, a crucial aspect of our design of the deformation is the adoption of a model in which contributions that may spoil the direction on a particular threshold singularity are excluded (i.e. "anti-selected"), as opposed to a model that enables (i.e. "selects") the correct contributions on the particular thresholds they are designed for.
We illustrate more specifically how an anti-selection model is preferable to a selection one by highlighting the shortcomings of the latter when applied to the previously introduced Box4E configuration whose four E-surfaces are shown in fig. 5 in the k z = 0 plane. The "selection" model would in this case amount to combine all four radial fields as follows (the discussion of the analogous construction of ref. [93] that involves normal fields would be -12 - where the selection function 2 simply is one minus the anti-selection function T (η bi ) defined as follows:T where M is an adjustable free parameter, and p 0 i −p 0 b is the length of the major axis of the Esurfaces η bi , which provides a measure for the size of the E-surfaces. Another possible choice is to substitute the normalisation , which is the minor axis length of the E-surface. The choice of M provides an estimate of how rapidly T (η) saturates to one when k is further away from the surface η bi .
The deformation of eq. (3.13) stemming from the selection model is problematic for mainly two reasons: • On the threshold E-surface η 12 , the deformation receives contributions mostly from v s 124 and v s 213 (which do satisfy the causal prescription) but also from v s 134 and v s 342 2 The selection function chosen in ref. [93] is an exponential Gaussian of adjustable width A bi :  Figure 5: A correct deformation direction with functional form described by eq. (3.17) for Box4E using four sources which are excluded on those E-surfaces whose interior does not contain the source. The right plot is a zoom-in on the central region.
(which may not satisfy the causal prescription) since the suppression factor induced by their respective selection function is small on this surface, but not zero. This implies the necessity of fine-tuning the suppression parameters which may be a difficult task when E-surfaces with very different causal constraints lie close to each other.
• On the intersection of two E-surfaces, for example ∂η 12 ∩ ∂η 13 , three of the four radial deformation fields v s 124 , v s 213 and v s 134 are active without any suppression, even though only v s 213 is guaranteed to be correct on this particular intersection.
One may think of alleviating the intersection problem by simply removing such intersections from the selector function applied to the deformation sources that are invalid: However, this solution is again not exact since even thoughT (η 42 )T (η 43 ) andT (η 43 )T (η 42 ) are small quantities on ∂η 12 ∩ ∂η 13 , they are not identically zero. In fact, it is impossible to build a continuous selection function that identically vanishes on a particular intersection of E-surfaces while at the same time being identically unity when evaluated anywhere on one of the intersecting E-surfaces but outside of the intersection. The above shows that if the the contour deformation is required to be correct (i.e. independently of its parameters), the radial deformation fields must be combined using an anti-selection paradigm that also avoids referring directly to intersections of E-surfaces, -14 -since one cannot continuously (anti-)select them. In the example of Box4E, we achieve this by constructing the final deformation vector κ as follows: which exactly satisfies the causal requirements for k on ∂η 12 and ∂η 12 ∩ ∂η 13 : In general, the minimal set of sources required for constructing a valid deformation with this anti-selection model is obtained by determining the maximal overlap structure of the E-surfaces, which we will formally define in sect. 3 12 , η 42 }}. After the maximal overlap structure has been determined, one has to construct source points in the interior of each overlap listed in the maximal overlap structure. Details about our strategy for choosing these particular points are given in sect. 6.1. Now that we have introduced and illustrated the key concepts underlying our construction of a valid deformation direction, we formalise it for an arbitrary number of loops and legs.

General solution to constraints on direction
In the absence of UV and IR non-integrable divergent behaviours, E-surfaces are the only singularities in the space of loop momenta that need to be regulated by a contour deformation. In sect. 3.1, we have shown that we have to construct a vector field pointing outwards on every E-surface. In this section we study this constraint in more detail. We remind the reader of the simplified notation identifying ( k 1 , . . . , k n ) with k that combines all coordinates of the n-loop integration space. E-surfaces are the boundary of convex, bounded volumes. We write the E-surface manifold as ∂η and its convex interior as ∂ − η, that is: The radial field k − s centred at point s has a strictly positive projection on any normal to the surface if and only if it is inside the surface itself: In general, given a set of E-surfaces F and a point in their interior: -15 -then k − s F will have positive projection on all normal vectors of E-surfaces in F and thus satisfies the causal prescription for all E-surfaces in F . We call s F the source of the set F . The aforementioned construction of the deformation field k − s F provides a systematic solution to the hard problem of constructing a deformation vector on the intersection of all E-surfaces in F , where many causal constraints need to be satisfied simultaneously.
In order to extend the applicability of the construction, we need to generalise it to more than one set of overlapping E-surfaces. Given the set of all existing E-surfaces E, we define the overlap structure (3.23) Thus O contains all possible sets of overlapping E-surfaces. One can immediately conclude that, if a set F is in O, then any subset F ⊆ F is in O.
Since a deformation vector k − s F is not guaranteed to satisfy the causal prescription on any point on an E-surface in E \ F , one has to identify the sets of overlaps F 1 , . . . , F N such that, among the radial fields k − s F 1 , . . . , k − s F N generated by such overlaps, there is at least one satisfying the correct causal direction on any point on an E-surface and, especially, on any intersection of them. Such a set with the least amount of elements is referred to as the maximal overlap structure O (max) and does not contain any set of E-surfaces that is a subset of another set in O: The set O (max) is the minimal set that ensures that one can build the final deformation without requiring special treatment for the intersections of E-surfaces (i.e. (anti-)selection thereof). Determining the maximal overlap structure is a challenging problem and is discussed in sect. 6.1.
In order to construct the deformation field for E, each element F ∈ O (max) is associated to a source s F whose corresponding radial deformation field k − s F is imposed to vanish on any E-surface not contained in F . This task is performed by a positive, bounded and smooth anti-selector function g F satisfying the following constraints In practice, we build g F ( k) from the same E-surface anti-selector building block T (η bi ), already introduced in eq. (3.15): which can be combined as follows to build g F ( k): Equipped with this anti-selection, we can now define a deformation field κ F valid for all E-surfaces in F (and their intersections) which does not contribute (i.e. it is exactly zero) to the deformation applied on any E-surface in E \ F : where the overlap function α F ( k) is, for now, any positive function which is non-zero on any E-surface contained in F . The construction of the final deformation can now be completed by adding together all vectors κ F where F ranges through at least all the elements of the maximal overlap set. We are now ready to write down a complete deformation field which satisfies the causal constraints stemming from all E-surface, independently of any deformation hyperparameter: The above minimal deformation field is what we used at one loop throughout this paper, including for producing the results presented in sect. 7. As we shall see in sect. 3.3.1, beyond one-loop it becomes necessary to consider additional deformation fields to accommodate particular continuity constraints of the integrand. We stress that supplementing the minimal deformation with additional causal fields can be performed without spoiling the causal properties of the individual terms because of the nature of the anti-selector functions. In fact, the sum κ F + κ F of two individually valid deformation vector fields κ F and κ F is also causally correct. More precisely, thanks to the anti-selection functions contained in κ F and κ F , we have that their sum is: • correct for k lying on an E-surface η in F or an E-surface η in F , but not on any intersection of η and η , that is on all points • exactly zero on the above-mentioned intersections as well as on any surface η not in F nor in F , that is on any point thus ensuring that κ F + κ F also satisfies all causal prescriptions if the deformation fields κ F and κ F already do. Another example of a deformation field that can be added is the sum of all appropriately anti-selected normal vectors of each E-surface. Thanks to this additive property of anti-selected deformation fields, one particular generalisation of eq. (3.29) is obtained by adding additional support sources from a set O of overlaps taken from the set The dependence of κ O ( k) on O underlines the aforementioned fact that adding to the minimal deformation vector -that is, the one constructed from O (max) -any deformation vector constructed from an extra overlap F ∈ O cannot spoil the causal constraints already satisfied by κ ∅ . More generally, it is also possible to add multiple radial fields generated by several sources from the same overlap F , although this is equivalent to adding a single radial field stemming from a different source in the same overlap. Adding support sources may improve numerical convergence and we intend to explore this possibility more systematically in future work. The particular strategy for selecting a near-optimal source point s F within a given overlap F is an implementation detail that we will discuss in sect. 6. The next section turns to the problem of assigning the correct normalisation to the deformation field constructed in this section. In particular, we will derive a necessary expression for the prefactors α F ( k).

General solution to constraints on magnitude
Once a procedure is established for constructing the correct deformation direction for a generic multi-loop integral, it remains to investigate conditions on the magnitude of this deformation. When writing the deformation vector field as λ κ( k), determining the normalisation of the deformation amounts to setting the value of λ. Constraints on the magnitude can be formulated locally for every k and can thus be satisfied by scaling parameters that are a continuous function of loop momenta λ = λ( k). For numerical stability it is typically advantageous to set the scaling parameter and the overlap function as large as possible while still satisfying the constraints.
The magnitude of the deformation is bounded by three conditions in the LTD framework: • the continuity constraint (sect. 3.3.1), • the expansion validity constraint (sect. 3.3.3), • the complex pole constraint (sect. 3.3.2).
Scaling parameters satisfying each of these constraints individually are denoted by λ cc ( k), λ e ( k) and λ p ( k) respectively. An overall scaling function λ( k) satisfying all three constraints can then be constructed as where λ max ∈ (0, ∞) is the maximum allowed value of the magnitude of the deformation. Although λ max is effectively a hyperparameter and thus subject to optimisation, the correctness of the deformation is independent of it. All the results presented in this work have been obtained by setting λ max = 10.
We will see that the continuity constraint also imposes conditions on the overlap function α F ( k) and the choice of overlap set O for eq. (3.32), thus arriving at the final expression for κ( k) that we will give in eq. (3.48). Our final expression of the contour deformation is then:

Continuity constraint
The request that the integrand is continuous on the contour adds constraints to the deformation vectors that have to be satisfied for all values of k, and specifically require that the argument of any square root appearing as energies of any on-shell particle never crosses the negative real axis, consistently with the choice of the principal square root branch. The energy can be written as a function of k − i κ: and thus the requirement of integrand continuity imposes that for any value of k and κ: Consider now a small ball centred at k * with q j ( k * ) = 0: then κ( k) has a constant direction throughout the infinitesimal volume of the ball (unless κ( k) ∝ k − k * ). Since q j ( k) spans all possible directions in this neighbourhood, it implies that there is always a continuous set of points containing k * and such that j . Instead of imposing this constraint on this continuous set of points only, we instead impose it everywhere, resulting in the following stronger (and simpler) version: which restricts the argument of the square root to lie in either the first or fourth complex quadrant. At one loop, given that q j ( k) = k + p j , this constraint can be satisfied by just using for the deformation from eq. (3.29) a scaling which imposes the deformation to always be lower in magnitude than E j ( k), ∀j, that is where cc is a parameter that we set to 0.95. The only problematic points are when a focal point of a massless internal propagator j, i.e. a solution of the equation q j ( k * ) = 0, j ∈ e, coincides with a point on another E-surface. According to eq. (3.38) this implies that E j ( k * ) = 0 and thus λ cc ( k * ) = 0, although the point is also located on an E-surface and thus requires a non-zero deformation. However, these points can be shown to be specific to the frame of reference initially chosen for the calculation and can be easily removed with a Lorentz boost (see section 5.1).
For multi-loop integrals satisfying the continuity constraint is not straightforward; indeed, consider an existing two loop surface equation for a massless diagram It admits as a solution the point ( a continuity constraint as in eq. (3.38) scales the deformation to zero, although the point itself is on a singular surface, and thus requires deformation. Strictly speaking, this dilemma is absent for diagrams with only massive internal propagators, as the masses act as regulators (i.e., Re[E 2 j ] > m 2 − Q j ( κ) 2 ) and forbid the deformation to be scaled to zero. However, in such cases a small mass imposes an unnecessarily strict constraint on the deformation in the neighbourhood of the corresponding focal point.
In order to remedy this problem, we observe that, given any proper subset of c ⊂ b of a loop momentum basis b, there is a proper subspace of the space of loop variables such that Q j ( κ) = 0 ∀j ∈ c, since the system is not full rank. This can be used to construct deformation vectors satisfying all causal constraints and branch cut constraints simultaneously on the portion of E-surfaces which lie on the subspaces q j ( k) = 0 ∀j ∈ c.
imposes conditions on s which make the radial field k− s automatically satisfy the continuity constraints in the neighbourhood of the subspace. The source determined this way is now partially constrained by the request that it satisfies the continuity condition without the use of a function directly suppressing the radial field on the subspace q j ( k) = 0, ∀j ∈ c.
One can now try to construct a deformation vector from sources satisfying eq. (3.40), by additionally imposing it has a causal direction on any E-surface when restricted to the subspace itself. More specifically, given the restriction of the E-surface to the subspace identified by c, the overlap structure is restricted to this subspace as well and can be defined as which is contained in the original overlap structure, that is O c ⊆ O. Given any element F ∈ O c , one can thus obtain a source s c F that satisfies the following convex constraints: Therefore, one can define a radial field k − s c F which will be non-zero on the subspace identified by c while still satisfying the continuity constraint and providing a causal direction on the portion of the E-surfaces in F and their intersections contained in the subspace identified by c. In order to not spoil causality outside the subset overlapping E-surfaces as contained in the subspace we will use a properly anti-selected deformation vector As before, κ c F will not violate causality constraints outside of the subspace, since the antiselector function g F will take care of setting the deformation to zero on E-surfaces corresponding to different overlaps in the subspaces characterised by c and all the E-surfaces not -20 -appearing in the subspace. Analogously to sect. 3.2, one can define the maximal overlap set in the subspace c (3.45) and thus construct a causal deformation vector when restricting integration to the subspace is exactly the deformation constructed in eq. (3.32) from the overlap structure obtained in the subspace identified by c, with all the overlap functions α c F ( k) chosen equal to a single function λ c ( k), which ensures that κ c ∅ ( k) satisfies the continuity constraint on any subspace different than c. That is: In order to construct the final multi-loop deformation vector field, it is necessary to associate a deformation vector to each strict subspace c ∈ P = b∈B P(b) \ {b} , where P(b) is the power set of the loop momentum basis b. We finally obtain where g F ( k) is the previously defined anti-selector function. Observe that eq. (3.48) is equal to eq. (3.29) at one loop since P = {∅}. Furthermore, since eq. (3.48) can be constructed from eq. (3.32) by setting it immediately follow that κ( k) is a causal deformation vector. One can observe that in the limit q j ( k) → 0 the deformation satisfies the continuity constraint Q j ( κ) 2 < q j ( k) 2 without necessarily being identically zero. We stress that, although the continuity constraint is satisfied on all subspaces and neighbouring points, there is no insurance that it is still the case away from it. Thus, as already mentioned, the final deformation vector must be given an overall scaling factor: which is now not suppressing the deformation to zero on subspaces. This concludes the construction of a general contour deformation which works both in the case of massive or massless propagators, satisfying all causal constraints.

Complex pole constraint
The analytically continued LTD integrand is singular at complex locations other than the real location of thresholds. These complex poles must not be included in the region of space between the deformed contour and the real hyperplane for the final result to be correct. This is consistent with the request that the integral on the contour matches the original one defined on R 3n .
The approximate complex pole location can easily be found when the square roots of E-surfaces are expanded up to second order in κ and the truncated expressions for the real part and imaginary part are set to zero: where the sum runs over all square roots expressing the energies appearing in the surface η (see eq. (2.7)), with the following coefficients Eq. (3.51) can be solved in the variable κ ∈ R 3n , for given k, which provides a parametrisation of the singular surface for the analytically continued integrand. Any point satisfying η( k) < 0 will admit no solution since the triangle inequality ensures that a i c i − b 2 i > 0, whereas points satisfying η( k) = 0 will have κ = 0 as a unique solution: the latter poles are the original E-surface boundary around which there is initially an intent to deform. Writing κ = κ n κ , we find that for η( k) > 0 there is a (3n − 2)-dimensional set of solutions which entirely lies on the hyperplanen κ · ∇η = 0 and which is radially symmetric with respect to the origin. This is illustrated for a two-dimensional example in fig. 6.
Whether a pole is included within the contour can be established according to the following guiding principle: given a parametrised deformation vector κ( k), the deformation contour will flatten out to become the original real space as the magnitude of the deformation κ( k) is sent to zero. Thus, if a pole is contained in the region between the contour and the real hyperplane for a given κ( k), κ can be scaled down such that the pole is exactly on the surface.
The request that the contour does not include any pole thus translates into a set of allowed values of κ for the deformation contour: κ is an allowed value if rescaling it so that κ → λ κ, there exists no value of λ ∈ (0, 1] such that a solution of eq. (3.51) is exactly on the contour. This immediately allows to state that, for given k, any value of κ satisfying is not allowed. Once the contour is explicitly parametrised as k − λi κ( k), the constraint on the allowed values of the deformation can be dynamically satisfied by using the treatment of ref. [88], which can be applied to any quadratic equation in the scaling parameter characterising the location of complex poles. Specifically, this treatment allows λ to take a large value whenever the imaginary part of the complex-valued surface is reasonably high in absolute value, as in these cases the deformation κ is far from the hypersurface orthogonal to the normal, which contains all the poles and forbidden areas. When κ approaches the surface orthogonal to the normal field, its value is constrained to yield a positive value for the real part of the surface η( k). In this way, the forbidden region eq. (3.51) is never reached. More specifically, given one has that there is no value of κ such that eq. (3.53) is satisfied if Finally, one can calculate and collect a scaling parameter λ η for each existing or nonexisting, pinched or non-pinched E-surface, and write It is important to include non-existing E-surfaces, as they may still have complex solutions.
It is particularly illuminating, in order to understand the relevance and location of the complex poles, to observe how the zeros of the original E-surface equation morph into the zeros of the real part of the complex valued E-surface equation. The location of the "displaced" threshold is implicitly determined through the equation This implicit equation defines a surface which is in general very different from the original E-surface, although it is clear that in the limit λ → 0, the two surface equations will be the same (see sect. 3.3.4 for visualisations). In the second order truncation in λ, it is also clear that the interior region of the displaced surface will necessarily contain the interior region of the original E-surface, since c i a i − b 2 i > 0, ∀ k. A rough bound on the volume of its interior region can be obtained by truncating the expansion of the square root to the next-to-leading order in the real part and requiring the correction to be smaller than (3.59) Figure 6: On the top left, an E-surface with its own normal field in ( k x , k y ) space. Three points, one in the interior of the E-surface (purple), one on the E-surface surface (blue) and one on the exterior of the E-surface (orange) are highlighted. In the other three pictures, one can find, for each of the highlighted points, the ( κ x , κ y ) space showing the forbidden line stemming from eq. (3.53) as well as the region allowed by the scaling of eq. (3.57) which guarantees that the deformation does not cross complex poles.
This equation can thus be used to provide an upper bound for the volume of the displaced threshold, in the form of another E-surface with the same focal points and larger constant term.
It is interesting to note that the real part of the complex-valued E-surface equation is negative in the interior region of the displaced threshold, and positive outside. It means that -24 -no forbidden values of the deformation can be crossed in the region outside the displaced threshold. However, inside the original E-surface, no pole is allowed. Thus, the region of loop momentum integration space which may lead to forbidden values of the deformation (when there is no appropriate dynamic scaling) is all contained between the original Esurface and the displaced threshold. An example of this behaviour is shown in fig. 9.

Expansion validity
The causal constraint on the direction and the complex pole constraint are formulated in the limit of a small deformation vector norm κ . In this limit, the imaginary part of η takes an especially simple form, as it prescribes that the projection of the deformation vector on the normal of η must always be positive. Likewise, the complex pole constraint admits an especially simple and elegant solution when η is expanded up to second order. This constraint also concerns the magnitude of the vector. Consider the energy where a j , b j , c j are defined as in eq. (3.52). Observe that the chosen stronger version of the continuity constraint eq. (3.37) already imposes that b j < a j and c j < a j . Thus a way to ensure the feasibility of the expansion is through the same mechanism which ensures that no branch cut is crossed. A more systematic approach to the constraints on the expansion, however, is to ensure that the argument of the square root is small in norm which leads to the condition This is effectively equivalent to requiring that the square root is expanded when its argument is contained within a disc of radius th . The overall expansion validity constraint can be satisfied by setting λ e equal to the minimal λ j for all energies E j : Another approach is to directly compare higher-order corrections to the leading order terms in the expansion. The odd orders are imaginary, whereas the even ones are real. The expansion to third order reads which, when compared with the expression -25 -yields the relation whose significance relies on the fact that suppressing the importance of the next-to-leading order with respect to the leading order of the expansion of the imaginary part also achieves the same for the real part. Suppression of this ratio can be obtained by imposing This shows that the choice λ 2 < a j 2b j makes the next to leading order contribution to the imaginary part dominate over the leading order when b j is small with respect to a j . As a consequence, the choice of the scaling of the deformation is constrained by the condition that The practical advantage of eq. (3.61) is that it is true to any order in the expansion, while its downside resides in a non-obvious interpretation of the expansion parameter . On the other hand, while eq. (3.67) only considers terms up to third order and does not account for the relevance of higher orders, it constrains the corrections to the imaginary and real parts simultaneously and consistently with only one expansion parameter. This parameter signifies the relative size of the higher-order correction with respect to the leading one.
The most conservative approach is to impose both constraints, but in practice we found good results by imposing eq. (3.68) only, which is what we used for producing the results presented in this work.

Visualisation of the contour deformation and its effects
In sect. 3.1 we constructed and visualised the deformation vector field for a one-loop configuration with four pairwise overlapping E-surfaces, called Box4E. In this section we will study the interplay between the contour deformation and the integrand in more detail.
First, we investigate the properties of the contour deformation k − i κ. Various aspects of the direction of the deformation vector κ were already discussed in sect. 3.1. In this section, we highlight details about the deformation magnitude κ , specifically, the impact of the three conditions it is subject to, as laid out in sect. 3.3. The magnitude κ can be studied at various stages in the construction of a deformation that will eventually satisfy all physical constraints. In fig. 7 we break down the construction of κ into four stages: (a) The deformation vector is subject to none of the constraints described in sect. 3.3 and the deformation magnitude is therefore determined alone by the superposition of all radial source fields.
(b) We impose the continuity constraint, introduced in sect. 3.3.1. It guarantees continuity of the integrand, since branch cuts of the square roots involved cannot be crossed thanks to this constraint. (d) The volume enclosed between the real hyper-plane and the contour deformation must not include any of the pole located at complex-values of the loop momenta. In order to guarantee this, we impose the complex pole constraint discussed in sect. 3.3.2. It again limits the magnitude of κ.
After these four steps, the deformation vector field κ is such that the integral is welldefined and yields the physically correct result. In fact, an E-surface η that has real solutions k ∈ (R 3 ) n of the equation η = 0 when the deformation is inactive ( κ = 0), has no more real solutions when the deformation is active. We therefore visualise the effect of the deformation on the E-surfaces. The deformed E-surface η defines two regions of interest: the zeros of its real part Re η, and the zeros of its imaginary part Im η. In fig. 8 we display the two regions of interest one-by-one for each of the four E-surfaces. With respect to the smooth elliptic surface described by η = 0, when the deformation is switched off, the regions Re η( k−i κ) = 0 and Im η( k − i κ) = 0 can be seen as a displacement of η into complex space. It is crucial here that these two regions do not intersect. If they did, i.e. the real and imaginary part of the E-surface equations were simultaneously zero, there exists a solution to the deformed E-surface equation η = 0, which cannot be allowed by our contour deformation, since κ satisfies the complex pole constraint. To showcase this exact scenario, we refer to the sideby-side comparison in fig. 9, where we used two deformation vector fields κ, a correct one and one that is not subject to the complex pole constraint. Its effect is subtle in this case, as it moves the real and imaginary solutions only marginally but essentially, as it renders the integral divergent without it. We take a more detailed look at the region between the four E-surfaces of the Box4E, as displayed in fig. 10. It contains the full deformation vector field κ and the regions of vanishing real or imaginary part of the deformed E-surface.
As a third aspect, we discuss how the deformation magnitude κ affects the integrand. The connection between magnitude and integrand becomes apparent when studying these quantities on a line segment in integration space. This line segment is displayed in fig. 11. We annotated 12 features, where one of them is a focal point and the remaining ones are zeros of either Im η = 0 or Re η = 0 of the deformed E-surface η. In fig. 12 we report the deformation magnitude κ along this line. We see that on the focal point the continuity constraint sets the deformation to zero (feature 1). At the other features the magnitude constraints lead to a non-smooth behaviour in the deformation vector field. In fig. 13 we study the integrand along the same line. We observe that on the focal point (feature 1) the integrand is singular. This is an integrable singularity and can be removed by using multi-channelling in the cut energies (see section 5.2).
Finally, in fig. 14 we show a density plot of the real and imaginary parts of the integrand, as well as the regions, where the real or imaginary parts of the deformed E-surfaces vanish. The enhancements in the real or imaginary part of the integrand are directly related to the zeros of the imaginary part of the deformed E-surfaces. These enhancements are expected when the deformation vanishes close to an E-surface.     Figure 13: The real (blue) and imaginary part (yellow) of the integrand multiplied by the Jacobian of the contour deformation along a line segment (see fig. 11) on a symmetric log y-axis: At annotated point 1, the line crosses a focal point. There, the integrand has an integrable singularity. Features 2 to 12 are crossings with the line and the points, where either the real or the imaginary part of the deformed E-surface vanishes. On these intersections the deformation vector field κ is non-smooth (see fig. 12), which induces discontinuities in the Jacobian of the contour deformation.

Subtraction
In the discussion so far, we considered integrals that do not have singularities for loop momenta of large magnitude (ultraviolet (UV) singularities) or soft and/or collinear to external legs (infrared (IR) singularities). For practical applications, such as computing amplitudes of physical processes, this will not be the case, as individual diagrams can contain both UV and IR divergences.
After transforming the integrand using LTD, non-integrable singularities manifest themselves as pinched (squeezed) E-surfaces. For the case of Feynman diagrams with massless internal propagators, this will happen when one or more of the massless external legs become on-shell. It is however still possible to numerically integrate such integrals, provided that the non-integrable singularities are regulated first. In general this is achieved by subtracting from the integrand an expression that contains the same pinched E-surface(s) and that approximates the original integral in the limit where the singular surface is approached. If these subtraction terms (also known as counterterms) are significantly simpler than the original integral, one can integrate them analytically in dimensional regularisation and add them back to the final expression in order to recover the original integral, including all its poles in the dimensional regulator. In this section we start by presenting a novel method to regulate divergent scalar integrals at one-loop without the introduction of propagators linear in the loop momentum featured in ref. [81]. We then discuss the introduction of counterterms for physical amplitudes [82] where only one term is introduced to remove all IR divergences. This regulated expression can then be integrated using LTD and the contour deformation discussed in sect. 3.
Note that in this section we refer to the external momenta as p i for ease of reading.

Divergent scalar integrals
We start by investigating scalar integrals subject to IR divergences at one-loop. In general, it is convenient to express counterterms in terms of the same building blocks as the original integrand, namely quadratic propagators. This allows to use the LTD formalism that has been introduced for the case of finite scalar integrals. At one-loop, we will show that we can always achieve such subtraction using a linear combination of triangles built by a subset of the original propagators and with coefficients expressed in terms of the kinematic invariants s ij . Since the counterterms involve only propagators already present in the original diagram, they do not introduce any new E-surfaces.

General one-loop massless scalar integral
Let us consider an n-point function with all the internal propagators massless and with external momenta p j with p 2 j = m 2 j . We first consider the case where only one leg i is massless (m i = 0). As a consequence, the corresponding scalar integrand will develop a collinear singularity when the loop four-momentum k becomes collinear to the corresponding -33 - . (4.1) In the expression above (where we consider the loop momentum to flow clockwise) we can see how the integrand factorises in the collinear limit. The integration of this counterterm can be performed as shown in ref. [81]. The variable x is a function of the loop momentum and is defined as follows: The expression on the l.h.s of eq. (4.1) can be written in an integral form as follows: (4.4) The coefficient c i (x) that multiplies the bubble propagators corresponds to the remaining hard propagators with the loop momentum evaluated in the collinear limit: The limit shown on the right-hand side of eq. (4.4) could be used to build an IR finite expression by subtracting it from I n , however such a counterterm introduces propagators that are linear in the loop momentum. Linear propagators yield singular surfaces that are not akin to E-surfaces, implying that the general construction of the contour deformation presented in sect. 3 cannot directly control the properties of the imaginary part of the loop momentum on them. We leave the investigation of solutions for accommodating linear propagators to future work and for now aim at casting the subtraction terms c i (x) in terms of propagators already present in the original divergent one-loop integral. We start by considering all possible triangles that factorise the same divergent bubble in the collinear limit. This condition fixes two of the three propagators of the triangle to be the ones that become singular in a specific collinear limit, whereas the third propagator -34 -can be chosen to be any of the other ones appearing in the original n-point integral. All such triangles are: with periodic conditions on the loop momenta labels. In the collinear limit, each element T (i, j) factorises one hard propagator t ij whose expression reads: Note that each squared momentum in the denominator of our coefficient functions is linear in x because p i is on-shell, resulting in only one simple pole in the variables x. In order to cancel the divergences of the n-point function we need to find a linear combination of T (i, j) with coefficients a ij (x) that satisfies: We can multiply both sides of this expression by the denominator of c i (x) which is equal to the product of all the possible t ij with i = j. We then obtain a polynomial of degree Since we have (n − 2) degrees of freedom and we insist that coefficients a ij (x) are free of poles in x, one needs to involve all terms T (i, j) in order to solve the equation above (assuming all the poles t ij (x) are distinct). In particular, an explicit solution can be found by using the roots of the inverse coefficients t −1 ij : resulting in coefficients that depend only on the external kinematics. This procedure does not work in the case of degenerate (raised) propagators. This can be resolved by considering a subsetJ i ⊂ J i which contains only one member of each degenerate subset of propagators with multiplicity ν j for j ∈J i . Moreover, we need to generalise eq. (4.7) in order to support the degeneracy of the involved propagators. In the collinear limit, the linear combination of the elements of this set gives the same singularities as the original integral, provided that: In this case we have |J i | parameters a ij to constrain a polynomial of degree n with (|J i |−1) distinct roots. It is then clear that the coefficients a ij take the same values as those given in equation (4.7). From this point onward, we will only consider one-loop scalar integrals with non-degenerate propagators.
We are now equipped with a method that removes single collinear singularities from integrals with one off-shell external momentum by writing a linear combination of the triangular elements T (i, j). When more than one external leg has a vanishing mass, we can apply the same procedure for each of them. In this case, we have to be careful when one of the triangles appears in more than one regularisation. For example, when two adjacent momenta are on-shell at the same time, one has T (i, i+1) = T (i+1, i−1). In this kinematic configuration the corresponding coefficients will be same: when Thus, one has to be careful when summing the regulator corresponding to each of the massless external legs in order to avoid double-counting. We can write one general subtraction term, referred to as CT n , that can be used for any combination of on/off-shell external momenta of a scalar one-loop n-point integral: where we introduced the coefficients β i used to avoid double counting. Their expression is where we make explicit use of the fact that whenever p i and p i+1 are on-shell at the same time the two coefficients a ij coincide. Because the constructed collinear counterterms do not depend on the parameter x, they completely remove the singularities from pinched E-surfaces, implying that they regulate both collinear and soft divergences. As a consequence, we have that the integral I n − CT n is finite for all loop momentum configurations. The original expression I n can be recovered by adding back the integrated counterterms. The integrated counterterm consists of n(n − 3) distinct one-loop scalar triangles that are straightforward to compute analytically for general external kinematics using dimensional regularisation. We leave to future work the investigation of the possible multi-loop generalisation of this construction of counterterms that do not involve any propagators that are linear in the loop momenta.

Explicit example of subtraction for a divergent one-loop scalar box
For the four-point box topology with massless propagators, there are four counterterms since the sum in eq. (4.8) over the coefficients a ij is empty. Only the β i are present and -36 -take the following expression: : otherwise , (4.10) where s ij = (p i + p j ) 2 . In the particular case where all external momenta are massless and on-shell (i.e. p 2 i = 0), the final expression of the counterterms reads: which coincides with the results presented in ref. [81], in which this same expression corresponds to the counterterm built for the subtraction of soft singularities (and the authors also concluded that the counterterm cancels all IR divergences in that particular case).
In other cases however, and especially beyond one-loop, the counterterms from ref. [81] introduce linear propagators of the form of eq. (4.1).

One-loop amplitudes
The first physical amplitude we study pertains to the production of photons from the scattering of a quark and an anti-quark. For brevity, the order of the photons is kept fixed during this discussion, as performing the integration over all permutations of the final states does not add any complications. The tree-level contribution for qq → (N − 2) V is defined as where all the fermions are assumed to be massless and the coefficients C 0 , T 0 depend on the vector boson considered as a final state. If only photons are considered as final states such coefficients are given by: These formulas can easily be extended to the electroweak bosons W ± and Z by substituting the photon polarisation vectors with generic ones / ε i →/ ε which also encode the information about the axial and vectorial part of the corresponding boson: -37 -with projectors defined as In order to obtain a more general expression we will use this new definition for the polarisation vectors. In the case of photons, all the P i s are proportional to the identity matrix.
In order to compute the one-loop QCD correction to eq. (4.12) one needs to consider all possible insertions of a gluon along the fermionic line. The IR structure of the relevant diagrams features one or two pinched collinear singularities if the gluon is attached to one or both the external fermion lines, respectively. In the latter case, the diagram also features a soft singularity.

Counterterms
If the photons are physically polarised, the only pinched divergences contributing to the IR sector involve a gluon connecting one of the propagators of the tree-level diagram with the external quarks. There are no singularities originating from two internal quarks and an external photon meeting at a vertex and becoming collinear, since the numerator vanishes: Since the pinched singularities originate uniquely from insertions of gluons connecting an external fermion to an internal fermion, the Ward identity can be used to regulate all the collinear and soft divergences with a general counterterm. However, it is necessary to fix a consistent choice of routing for the loop momentum in order for cancelling divergences to be localised in the same region in momentum space, even though they belong to different diagrams. The general counterterm I IR reads: where This integration can be performed analytically using Feynman parametrisation, and we obtain: (4.20) Although subtracting eq. (4.17) from the original integrand allows to completely regulate IR singularities, the subtracted integrand is still divergent in the UV sector. This divergence can manifest itself locally, in spite of the integral itself being finite, either due to symmetries of the integrated expression or because the IR and UV poles cancel for integrals that are scaleless in dimensional regularisation. The behaviour for large momenta is inferred by the scaling of the integrand in these regions, and as a result all log-divergent triangles (one gluon, two fermions) and linearly divergent bubbles (one gluon, one fermion) that appear in the amplitude have to be regulated. The construction of the counterterm is done by taking the UV limit of each diagram by replacing where the only relevant momentum is now the loop momentum carried by the exchanged gluon. The bubble diagram has a leading UV divergence that is linear in the loop momentum. In the context of an analytic integration such contribution integrates to zero because of radial symmetry, although the integrand is locally divergent. It is therefore necessary to also regulate this leading UV divergence together with the sub-leading one obtained by computing the second order in the Taylor expansion around the UV approximation given by eq. (4.21). An explicit example of this subtraction can be found in appendix B, where eq. (B.11) represents the UV counterterm of a triangle and eq. (B.10) represents the counterterm of a bubble. The IR counterterm that we introduced is UV divergent and requires regulation as well. Its divergence can be expressed as as a triangle integral and can be subtracted by means of eq. (4.21). The combination of counterterms can be used to build a finite amplitude expression that can be integrated using LTD: The counterterm can be integrated analytically with the use of dimensional regularisation. In the UV contribution to the integrated counterterm we notice that the bubble and the triangle lead to the same value in norm and opposite in sign if constructed according to the substitution rule (4.21). Thus, the only remaining contribution is UV div.
The complete expression I CT can then be expanded in up to finite terms and be used to recover the original amplitude once combined with the value coming from numerical integration. The integrated counterterm for qq to photons takes the simple form: 25) where ln µ = log µ 2 −s 12 . Any dependence on µ UV has dropped from this final expression. As a consequence, the integration of the finite amplitude will also not depend on the choice of µ UV . This condition can be used as a further check for the proper cancellation of the divergences.

Ultraviolet behaviour
When integrating the LTD expression, one has to take into account that the superficial degree of UV divergence of each dual integrand is higher than that of the sum of its cuts. This is because once the LTD on-shell cuts of the residues are applied, every quadratic propagator scales as 1/| k| in the UV instead of 1/k 2 . As a consequence, contrary to the Minkowskian case, the addition of more fermion propagators to the diagram is not suppressing the scaling of the deformation in the UV sector: compared to the original scaling of the 4D integrand being Summing over all the different cuts will however recover the original scaling of k 2−N . If the dual integrand scales faster than 1/| k| in the UV, the numerical cancellation of large numbers becomes prone to numerical instabilities. One way avoid such numerical instabilities in the UV region is to approximate the integrand with a better behaved function in the corresponding sector, obtained by taking a UV approximation of the integrand. The most convenient choice is to replace all the propagators with a common UV one: (4. 28) This ensures that the approximating function only features a single dual integrand, which directly scales as the 4-dimensional integrand. The numerator can be left unchanged for this approximation. In section 4.2.3 we discuss the effects of this UV approximation.
-40 - The UV counterterms can be constructed as shown in sect. 4.2.1 for most integrals, but in the case of a bubble integral, the subleading logarithmic divergence must also be regulated. The relevant part of the approximation is shown below: (4.29) Since the UV counterterms have higher-order poles, the LTD formula shown in sect. 2 cannot be applied directly. We discuss how to apply LTD to integrals featuring raised propagators in appendix C.

One-loop amplitude for
We now study the specific case of the one-loop dd → γ 1 γ 2 γ 3 amplitude. The tree-level diagram of this amplitude is where the coefficients are given by (4.30) Fig. 15 shows all the diagrams involved in the one-loop QCD correction. Diagrams D1 -D3 and D7 -D8 are IR divergent: D1 and D7 are divergent when k is collinear to p 1 and D2 and D8 are divergent k is collinear to p 2 , whereas the diagram D3 is divergent in both cases and also has a soft divergence.
Despite the fact that the integrated amplitude is UV finite, the local behaviour of the integrand in the UV region needs to be regulated. This can be done by writing the corresponding counterterms for all UV divergent integrals, specifically D4 -D8.
In order to ensure that the cancellation occurring across diagrams at the integrated level are also reflected at the local integrand level for the whole amplitude, one must carefully choose the the loop momentum routing of each diagram so as to localise cancelling divergences in the same region of momentum space. The case at hand is quite easy in that regard, as one can choose the gluon line to have momentum k with momentum flow against the fermionic line for all the diagrams. Fig. 15 shows the different behaviours when approaching the soft, collinear, and UV limits. The different limits are approached by rescaling the loop momentum k by a factor δ for the soft and UV limit, while for the collinear limit we use the Sudakov parametrisation Figure 15: Diagrams contributing to one-loop QCD correction to qq → 3 γ amplitude.
of eq. (4.2) with y and k ⊥ rescaled by δ and √ δ respectively. The different asymptotic scaling δ 1 , δ 1 2 and δ −1 , prove that the divergences are properly subtracted. Despite the use of quadruple precision (f128) to rescue some unstable evaluation of the UV region, we see that the cancellations between dual integrands are broken around δ > 10 8 due to numerical instabilities. In fig. 17 we show how these instabilities spoil the final result in the case of double precision (f64) with and without the use of the approximating function discussed in sect. 4.2.2. In the latter case it is possible to push the instability in the far UV and reproduce the behaviour of the quadruple precision evaluation. Where the transition between the approximated function and the all-order amplitude expression occurs, one has has to ensure that the deformation goes to zero, since this region is not analytic. In both fig. 16 and fig. 17 the rescaled loop momentum is taken to be real and of the same order as s 12 .
-42 -    Figure 17: Behaviour of numerical instability in the UV due to imprecise cancellations between large numbers from each each dual integrand. The loop momentum k is rescaled by a factor δ and the real and imaginary part of the amplitude are presented with different precision (double and quadruple) and by expanding the expression around the UV limit as an approximation (see sect. 4.2.2).

Optimisation
In this section we present various optimisations that we have developed to improve the convergence of our numerical framework.

Lorentz invariance
The following two subsections are aimed at showcasing the wide range of simplifications made possible by leveraging Lorentz symmetry. Specifically, Lorentz symmetry can be used to both drastically simplify the E-surface overlap structure and eliminate fictitious accidental pinched configurations that may appear for specific external kinematics as a result of competing constraints on the deformation. Contrary to symmetry under the (spatial) SO(3) subgroup of the Lorentz group, invariance under boosts is not manifest in the LTD framework. Indeed, Lorentz boosts cause significant changes in the singular structure of the integrand and result in E-surfaces being rescaled and shifted relative to each other: the major axis length of an E-surface, being a linear combination of the energies of the external particles, is not a Lorentz invariant, nor is the distance between any pair of focal points, being a linear combination of the three momentum of the external particles. Conversely, some quantities are Lorentz invariant in the LTD framework: the number of E-surfaces, their existence condition, and some specific features of the overlap structure including, for example, the property of two E-surfaces sharing a focal point.

Simplified deformation contour for 2-point multi-loop integrals
A first use-case of the implict realisation of Lorentz invariance in LTD is found in the construction of a surprisingly simple integration contour applicable to any two-point function.
Since the original integral is Lorentz invariant, the single independent external momentum of a two point function can always be boosted in its rest frame. It follows that the spatial momentum shifts in all propagator momenta read where we recall that e identifies the list of edges of the loop graph. Equivalently, we can write q i ( k) = Q i ( k). A Lorentz boost thus allows to decouple components of k from the spatial part of the external momentum. This feature allows for a simpler deformation, characterised by the parameter λ ∈ (0, 1), as This deformation casts squared energies in a particularly simple form, from which follows that because λ < 1, the stronger continuity constraint eq. (3.37) is always satisfied, since the real part of eq. (5.3) is positive and that all focal points coincide with the origin thanks to eq. (5.1). And because λ > 0, the imaginary part of eq. (5.3) is positive as well. It follows that the causal constraints, imposed by LTD, are satisfied everywhere (except at the origin where the deformation scales to zero), since the deformation is guaranteed to never reach the forbidden areas presented in eq. (3.53). Therefore, the simple deformation vector field κ = iλ k with λ ∈ (0, 1), is correct for any two-point function, independently of the number of loops and internal masses.
We tested this deformation on a six-loop two-point ladder integral with two sets of kinematic configurations given by p 2 = 1 and masses m 2 j = 0, ∀j ∈ e, called K, and p 2 = 1, m 2 j = 0.1 ∀j ∈ e called K . We compared the m 2 = 0 numerical result against its analytical counterpart and verified that the procedure is correct. The results are reported in the following The same technique of adding a small imaginary part to the components of the loop momenta corresponding to zero components of all the external momenta can also be considered for the three-(four-)point function. However, in these cases there are only two(one) component(s) that can be set to zero through a boost. The possibility of integrating easily along loop momentum dimensions by adding a small imaginary part to a subset of the components of the loop momenta is the manifestation of a property of two, three and four-point functions already noted in ref. [97].

Example of overlap structure simplification for a 3-point 2-loop integral
In general, Lorentz boosts can be used to greatly simplify the overlap structure. For example, we find that the 1 → 2 kinematics of a two-loop ladder diagram with massless propagators (considered here for simplicity), can be written in the following form when boosted in the rest frame of the p 2 + p 3 system: with momentum conservation conditions yet to be applied to the energy components. Since in this case any E-surface features at most one focal point with a non-vanishing affine term p j , the origin k i = 0 lies within all E-surfaces. Indeed, all E-surfaces of this particular loop integral considered are which are all negative when evaluated at k i = 0, indicating that the origin is indeed in the interior of all exisitng E-surfaces. Similar arguments can be used to show that in a physical 2 → 2 process featuring n existing E-surfaces, at least n − 1 of them must allow for a point in the interior of all of them. The boost parameters can themselves be viewed as hyperparameters subject to optimisation and although it is beneficial to boost 2 → N kinematics in the rest frame of the collision, a systematic procedure that maximally optimises the choice of Lorentz frame is still missing.

Pseudo-pinches
Pseudo-pinches are singular surfaces at which competing causal or continuity constraints impose the deformation to be zero, although these configurations are non-existent in another frame of reference. They can be classified as follows: 1. Singular subspaces with |c| fixed loop variables and n − |c| unconstrained loop variables. When all loop momenta configurations k (c) satisfying the subspace constraints of eq. (5.6) happen to also lie on one particular E-surface η (so η( k (c) ) = 0), then no deformation will be allowed on that surface because of the continuity constraint of eq. (3.3.1). This situation is accidental as it only happens for particular kinematic configurations and, more importantly, for a particular choice of Lorentz frame. At one loop, this situation corresponds to a focal point being located exactly on an E-surface.
2. Intersections of two or more E-surfaces η 1 , . . . , η n at a point k such that ∃η i with ∇η i = − j =i α j ∇η j and α j ≥ 0. This typically happens when two E-surfaces are tangent. We stress here again that, in general, the normal ∇η i to an E-surface η i is a (3n)-dimensional vector.
We now illustrate these two different types of accidental pseudo-pinches at one loop.
Case 1 Let a focus be located exactly on an E-surface. Imposing that the contour does not cross branch cuts of on-shell energies of massless internal particles (using our stronger version of the continuity constraint), -46 -at the point q i ( k * ) = k * + p j = 0 implies that κ * = 0. However, since the point is located on a singular E-surface, k * ∈ ∂η, a non-zero deformation is required. In this case, the continuity constraint conflicts with the causal constraint. It can be argued that our continuity constraint is stronger than what is minimally required, but even weaker implementations must impose that κ * = 0 in some region containing the focal point.
Case 2 Now let two E-surfaces be tangent. Then, two causal constraints conflict at a point: the normal vectors to the two E-surfaces at the tangent points are opposite in direction, and thus no vector exists having strictly positive projection on both of them.
Both cases are problematic from a conceptual point of view, because they can correspond to kinematic configurations where the deformation breaks down. However, as mentioned earlier, the existence of these cases is accidental and specific to the chosen reference frame for the external kinematics. In both cases, there is an infinite number of infinitesimal Lorentz boosts such that in the boosted kinematics no focal point coincides with any E-surface and no two E-surface are tangent. This is especially clear in the case of causally connected focal points. In order to understand this notion, one can turn to the one-loop example of an E-surface η on which lies a focal point f (necessarily, the focal point f cannot coincide with one of the focal points of η). Now let q i ( k f ) = k f + p i = 0 be the equation defining the focal point and let f be a focal point of the E-surface satisfying the equation q j ( k f ) = k f + p j = 0. Now consider a boost sending the four-momentum p i − p j in its rest frame so that its only non-zero component is the time component. Obviously, this can only be done if p i − p j is timelike in which case the two focal points correspond to four-dimensional spacetime coordinates that are causally connected. In this frame of reference, the focal points f and f overlap and thus f can no longer be located on the surface of the ellipsoid, thereby avoiding the accidental pseudo pinch situation.
Similarly, consider two tangent E-surfaces, and choose one focal point for each Esurface, denoted by f and f , such that their distance in four-dimensional spacetime is timelike. It is now always possible to choose a frame of reference in which the distance between the focal points is zero. In this frame the two E-surfaces share a focal point and thus cannot be tangent.

Multi-channelling
Improving the numerical efficiency of the numerical integration amounts to finding techniques for reducing the variance of the integrand. Sharp local enhancements of the integrand, and especially integrable singularities, induce a large variance and can significantly deteriorate the numerical integration. At best, such peaks make the Monte Carlo (MC) integration converge slowly and at worst they yield an unstable central value, as well as an unreliable estimate of the MC error.
In general, adaptive importance sampling can adjust well to integrands with large variances, provided that their enhancement structure aligns with the integration variables. However, when the Monte Carlo integrator underestimates the variance of the integrand in -47 -some regions of the integration space during the first iterations, it can incorrectly neglect these regions in further iterations. In such cases, the estimate of the integral will be unreliable, even though the error suggests otherwise. Even though increasing the number of sampling points in the first iterations can help mitigate this problem, it slows down the integration and reduces the predictive power of the numerical integration. It is therefore best to first pre-process the integrand so as to remove its sharp enhancements, which is possible when their location and approximate functional form is known. In this section, we show how this improvement can be systematically implemented for the LTD expression, using a technique known as multi-channeling which is commonly used for improving numerical integration in various contexts.
We can write the integrand stemming from the n-loop LTD expression as where each dual integrand Res b [f ] features sharp peaks resulting from each propagator put on-shell. Each of these peaks is an integrable singularity when the corresponding propagator is massless. These enhancements for each residue have the following functional form: In order to take advantage of dual cancellations, i.e. the local cancellations of singularities on H-surfaces among summands of the LTD expression, the dual integrands have to be integrated together using a unique parameterisation. We must therefore consider the complete integrand which features the following peak structure In a multi-channeling approach, we seek to flatten these enhancements by first inserting the following expression of unity in the integrand: and then splitting up the sum in the numerator into |B| channels, thereby defining an integrand for each channel identified by a basis (or equivalently spanning tree) b ∈ B, whose expression reads: We observe that each channel still features peaks, but only those specitic to b. This opens the possibility of choosing a different parametrisation for each channel, selected so that its Jacobian flattens its enhancement i∈b E −1 i . We note that a similar multi-channeling approach was used in refs. [72,75]. Thanks to the continuity constraint discussed in sect. 3.3.1, the denominator of the multi-channelling factor does not introduce new integrable singularities when computed with our choice of contour deformation. More specifically, the integration measure from the spherical parametrisation of the loop momenta in the basis b reads 3 : where we introduced the shorthand notation We can now choose to integrate each channel C b separately 4 and use for each the specific parametrisation of eq. (5.13). At one loop, these different parametrisations only differ by a shift of the origin whereas beyond one loop, they also amount to a change of basis in which the loop momenta are expressed. The resulting integral for each channel then reads: where each of the two factors building the integrand is now free from integrable singularities (or strong enhancement in the case of massive propagators) coming for the cut propagator. The original integral is then computed as the sum of |B| channels The effects of multi-channeling are shown in fig. 18, where the peak due to the crossing a focal point is removed.
We note that this multi-channeling approach can be further developed by considering additional channels related to other enhancements coming from E-surfaces and/or infrared limits for example. We leave this investigation to future work.

Numerical implementation
In this section we discuss various details of our numerical implementation, such as the most challenging aspects associated to the construction of the deformation contour, the evaluation of the Jacobian and consistency checks that are essential for verifying the correctness of the integration contour and guaranteeing the stability of the evaluation of the integrand. Since the triangle has three dual integrands, the LTD integrand I has three integrable singularites, one for each energy E i = 0, i ∈ e. For both integrands, the singularity at E 3 = 0, i.e. when k = 0, vanishes when parameterised in spherical coordinates centered at k = 0 because of the integration measure. The line along − p 1 goes directly through the singularity at E 1 = 0, i.e. when k = − p 1 and past the one at E 2 = 0 (small bump only since the direction used for this plot is not p 2 but p 1 ) of the LTD integrand. In the channel C {3} these two enhancements are flattened and become non-vanishing constants thanks to the multi-channel factor. We observe that at k = − p 1 the channel is not differentiable (as well as at k = − p 2 ).

Source determination
Determining the maximal overlap structure requires testing whether there is a point in the interior of a given set of E-surfaces. This problem is convex and, in particular, can be written as a second-order cone program (SOCP). We have used the convex constraint problem rewriter cvxpy [98] with the ecos solver [99] as a backend to construct a program that ascertains whether a given set of E-surfaces overlap. Given the aforementioned program, determining the maximal overlap structure O (max) of eq. (3.24) is still an NP-hard problem, as the set of possible overlap configurations is exponential in the number of E-surfaces and any algorithm devoted to the determination of O (max) will have a worse-case complexity that renders it prohibitively slow. In practice however, the class of problems of interest generally features a limited amount of overlapping regions which are shared by many E-surfaces. Indeed, many E-surfaces share one or more focal points, and thus naturally have the focus as a shared interior point. As a consequence of these facts, the algorithm should be constructed so as to take advantage of this heuristic my exploring solution in a top-down order; that is starting with the assumption that all -50 -E-surfaces overlap. If all E-surfaces are not in one overlapping set, one E-surface is removed in all possible ways and the test is performed again. Once an overlap is found involving N particular E-surfaces, then the 2 N − 1 subsets of this set never need to be tested again. In order to prevent a combinatorial blow-up, a list of all possible pair-wise intersecting E-surfaces is constructed and used to filter many options when constructing viable subsets. This additional improvement to the heuristic was key in rendering our implementation fast enough for problems with more than 30 E-surfaces, as generating all 2 30 options is too slow. In practice, the refined algorithm takes only a few seconds to find the solution in the majority of cases. It therefore yields negligible overhead in comparison to time spent in the numerical integration. We note however that for cases involving or more that 40 E-surfaces, it may happen that when our heuristics are not well satisfied, our algorithm cannot determine the maximal overlap structure within any reasonable amount of time, as it happened in the case of the loop integral 7.2L8P.K1 * for which we could then not show results.
Once the maximal overlap structure is determined, one must find a point inside each overlap with the extra property to be optimal from a numerical convergence point of view. This optimality condition can loosely be approximated by requiring the point to be as far as as possible from all the E-surface defining and enclosing the overlapping volume. The resulting set of point constructed in this manner will serve as the set of deformation sources. The furthest away a source s is from all surfaces in the overlap set, the less tangential the deformation k − s will be when evaluated on the surfaces themselves. For higher-loop cases, the source location is possibly subject to extra requirements due to the continuity constraints within a particular subspace given in eq. (3.40).
To approximate the optimal centre of the overlap region, which is related to the Chebyschev centre of a convex region, one can solve the convex constrained optimisation problem of maximising the radius r under the constraints that the points s ± rê (j) i all lie inside all E-surfaces η ∈ F for every Cartesian direction in 3n dimensions, {{ê ≤ 0, ∀i = 1, 2, 3 ∀j = 1, . . . , n ∀η ∈ F Imposing the extra subspace constraints of eq. (3.40) is most conveniently done by performing a basis change. For example, for given linear constraints k 1 = p 1 and k 1 + k 2 = p 2 on vectors ( k 1 , k 2 , k 3 ), the following system of equations allows to identify the subspace satisfying the constraints and its orthogonal complement where ker(C) is the kernel of the constraints C, (0, 0, 1) in this example. The inverse of the system presented above allows to rewrite the E-surfaces in terms of fixed momenta p 1 , p 2 and the source variable s 1 . In this particular subspace example, there remains only three degrees of freedom for setting the source, so that only three canonical directions e i , whereas the original centre finding problem cast without change of basis would require all nine (3n).

Parameterisation
The numerical integrator Cuba [100] that we use to produce our results generates points in the unit hypercube [0, 1] 3n . These points have to be transformed to R 3n where they then correspond to a particular real-valued sample configuration for the spatial part of the the loop momenta. Our code provides options for Cartesian maps and spherical maps with hyperbolic and logarithmic scaling for the conformal mapping from [0, 1] to (−∞, ∞). For the results in this paper we used the following spherical and hyperbolic transformation that map each triplet of input variables (u 1 , u 2 , u 3 ) ∈ [0, 1] 3 to a configuration of the spatial part of one loop momentum k: where E cm is the centre-of-mass energy of the decay or scattering kinemtics, and b is a scaling parameter that regulates how much the integrator probes the ultraviolet region. Our default value for b is 1.

Deformation Jacobian
The contour deformation k → k − iλ( k) κ is effectively parametrised by the real part of the loop-momenta. Determining the resuling Jacobian of this parametrisation analytically is difficult due to off-diagonal contributions in the Jacobian matrix from the generally complicated analytical expression of the deformation magnitude λ( k). In order to bypass this inconvenience, the exact Jacobian is calculated numerically using automatic differentiation. This technique is commonly used in machine learning algorithms, such as neural networks. Performing the computation with dual numbers where the dual components i are subject to the truncation rule i j = 0, yields the partial derivatives ∂k ∂k jo as the coefficient of k jo . In our Rust implementation, all routines are generic over floating-point-like types (such as a double-precision floating point number). Since a dual number behaves like a floating point number, the promotion of the arithmetics to dual number can be done transparently from the perspective of our core routines implementing the LTD logic.

Consistency checks
In order to assess the numerical stability of each evaluation, each Monte Carlo sample point is evaluated on numerically different but analytically equivalent integrands, taking advantage of the manifest invariance of the integrand under rotation of the spatial part of every momentum involved (for example, the external momenta, the loop momenta and the sources). If the evaluation of the LTD integrand of a spatially rotated configuration significantly differs (in terms of a sensible adimensional threshold) from the original one, the point is deemed unstable, and we attempt to rescue it by repeating the same exact procedure in quadruple precision. If an unstable point is then considered stable in quadruple precision by performing the same test, then the quadruple-precision result is returned to the integrator. Instead, if the point is still deemed numerically unstable, we set its weight to zero. In practice, even for the more challenging integrals, less than one sample point in a million is numerically unstable in quadruple precision. Furthermore, these exceptional unstable points are often deep in the ultraviolet region and evaluate to values far below the result of the integration and they can therefore safely be set to zero. We note however that the implementation of a quadruple precision rescuing system was necessary for obtaining many of the results presented in this publication, especially for the computation of amplitude where the ultraviolet behaviour is more relevant (see sect. 4.2.2).
The correctness of the complex contour deformation is verified by sampling random points on E-surfaces and ensuring that the causality constraint is satisfied. Since finding a parametrisation for E-surfaces is difficult at higher loops, it is more effective to use a bisection strategy to sample points on the E-surfaces. The bisection strategy must be seeded by one point inside the E-surface and one outside. As E-surfaces are bounded, finding a point in the exterior of them is trivial and the most straightforward choice of point in the interior is any of the two focal points of the E-surface. The convexity of E-surfaces then ensures that a unique (correct) solution will be found by the bisection algorithm and that all points of a given E-surface can be reached by our approach simply by varying the choice of exterior point.
To verify the validity of the LTD expression, the occurrence of dual cancellations is explicitly verified. A similar bisection strategy is used to find a point on an H-surface. Then along the bisection line, the LTD integrand is evaluated on points iteratively closer to the H-surface. If the slope of the interpolation between these points is below a chosen adimensional threshold, the dual cancellation is considered successful. The same setup is also used to verify if the local counterterms used to subtract IR-divergences have the correct scaling behaviour (see section 4).

Results
The aim of our work is to provide a numerical loop integration technique based on Loop Tree duality which is both robust and generically applicable. It is therefore crucial to accompany the formal derivation of a valid deformation carried out in sect. 3 with illustrative applications that can demonstrate the correctness of the numerical method as well as its practical efficiency. This will be explored in sect. 7.1. We present our numerical results -53 -obtained when applying our LTD formulation together with local subtraction counterterms to compute one-loop scalar topologies in sect. 7.2 and to compute amplitudes for the ordered production of photons from a fermion line in sect. 7.3 .

Multi-loop finite integrals
To demonstrate the practical efficiency and correctness of the deformation, we explore in tables 1-8 a variety of kinematic configurations and many different scalar integral topologies featuring up to four loops (and up to six for cases not necessitating a contour deformation) 5 , yielding different combinations of number N E of unique singular threshold E-surfaces and number N S of necessary deformation sources. We also indicate the number of dual integrands in the LTD expression of eq. (2.5) in the column labelled N c ; it corresponds to the number of spanning trees of the topology and also to the number of integration channel it would feature when adopting the multi-channeling procedure discussed in sect. 5.2 (which we do not use in this section, unless otherwise stated).
We also report a shortened representation of the maximal overlap structure O (max) as a list L max where each entry corresponds to the number of E-surfaces contributing to each maximally overlapping set F contained in O (max) . We report the discrepancy of our numerical LTD result w.r.t the reference value, relative to each other (∆[%]) and relative to the Monte-Carlo error (∆[σ]) reported by the implementation in Cuba [100] of the Vegas [101] integrator 6 . Unless otherwise stated, we consider different fixed statistics of 3 · 10 9 , 1 · 10 9 and 0.5 · 10 9 Monte-Carlo sample points for each of the one-, two-, three-and higher-loop integrals computed 7 . For some of the one-loop results (e.g. 1.1L5P.V and 1.1L6P.IX), the real part is accidentally small compared to the imaginary part and since the variance of the LTD integrand is of the same order for both phases, we find it relevant to also indicate in the last column of the results table the relative discrepancy of our LTD numerical result on the modulus of the complex-valued benchmark result (∆[%]| · |). The timing per PS point t/p is reported in microseconds, as measured on a single core of an Intel Xeon CPU E5-2650 v4 @ 2.20GHz CPU. Throughout this section and unless otherwise mentioned, we keep the deformation hyperparameters fixed to their default values of th = 0.3 and M = 0.07. These defaults are typically different from what would be the values optimised for each kinematic configuration and/or topology tested, but in this exploratory work we refrained from systematically fine-tuning hyperparameters so as to prevent any bias in our results and be able to fairly showcase the robustness of our approach. However, we will later show two examples where the results from specific integrals could be significantly improved by adjusting the value of the hyperparameter M . Finally, the reference result for all one-loop integrals presented in this section, as well as for the one-loop amplitude computed in sect. 7.3, is obtained from the One-Loop Provider MadLoop [102,103]. MadLoop uses 5 The exhaustive details (incl. kinematics) necessary for reproducing the results of each integral presented in this section is given in the ancillary material. The integral normalisation matches that of eq. (2.1). 6 Similarly to the findings of ref. [80], we also find significantly more accurate and precise results using the Cuhre integrator at one-loop. The results with this integrator are however significantly worse beyond one-loop. For the sake of simplifying the comparison of our results across loop counts, we only report results obtained with the Vegas integrator. 7 With typically n_start ∼ 1% of n_max of and n_increase ∼ 0.1 % of n_max in Vegas.
In table 1, we present results for one-loop five-and six-point scalar integrals for handcrafted kinematic configurations that correspond to many qualitatively different maximal overlapping situations. We also include the result for the four-point one-loop integral 1.Box4E which we used as an example throughout this work. The relatively good sub per-mil accuracy obtained for this integral may be surprising in regard to the complexity of the corresponding LTD integrand, depicted in figs. 13 and 14. Comparing the Monte-Carlo accuracy and precision obtained for all integrals of table 1, we observe the general trend that the convergence mildly degrades with an increase in the number of deformation sources and the number of unique threshold E-surfaces. However, the dominant factor appears to be the shape of the threshold surfaces, which become more elongated as the masses of the external momenta decreases or, more in general, when the hierarchy between the relevant scales in the scattering considered becomes more pronounced. The integrals 1.1L6P.VII and 1.1L6P.VIII are a prime example of this observation as the Monte-Carlo accuracy of the latter integral is much worse despite featuring the same number of unique E-surfaces and deformation sources as the former. Indeed, the external kinematics of integral 1.1L6P.VIII yield E-surfaces of very elongated shapes, as hinted by the corresponding maximal overlap structure L max = [3,5,6,7] where one deformation source involves only three out of the total of ten unique threshold E-surfaces. Fig. 19 shows a rendering of the E-surfaces from both integrals 1.1L6P.VII and 1.1L6P.VIII, which clearly highlights their differences in shape and maximally overlapping regions.  Figure 19: The singular E-surfaces from the two six-point one-loop integrals 1.1L6P.VII (left) and 1.1L6P.VIII (right) with different kinematics yielding drastically different maximally overlapping regions. In both cases our construction of the deformation is generated from the combination of four radial fields with sources indicated by blue dots. Additional support sources can potential improve on the worse convergence of integral 1.1L6P.VIII. . On the other hand, the column L max in the multi-loop case still refers to the cardinality of the sets in O (max) (that is, the maximal overlap structure obtained in the absence of any focal point constraints). Furthermore, beyond on loop, the number of channels (i.e. number of dual LTD integrands) N c is no longer equal to the number of propagators, but instead corresponds to the number of spanning trees which is a quantity specific to each integral topology. Integrals 2.2L6P.a.I to 2.2L6P.f.I reproduce results from ref. [89], in which the authors perform a direct integration in four-dimensional Minkowski momentum space. We investigate the exact same decay kinematic configurations as the ones considered in that work, which are numerically well-behaved and yield results that are pure phases. We also obtained independent reference results for these two-loop six-point integrals using an alternative numerical computation using pySecDec [62] and we find only small tensions between all three results.
The multi-loop ladder four-point integrals (2.2L4P.c.I, 3.3L4P.I, 3.4L4P.b.I, 3.5L4P.I and 3.6L4P.a.I) are known analytically for massless internal lines [96], and a generalisation to M xN fishnet topologies (of which integrals 3.4L4P.a.I and 3.6L4P.b.I are two examples) was recently carried out in ref. [110]. We stress that the five-and six-loop integrals 3.5L4P.I, 3.6L4P.a.I and 3.6L4P.b.I are computed for external kinematics yielding no threshold singularities such that the integration can be be performed without any contour deformation. Furthermore, for these integrals, we used the multi-channeling treatment discussed in sect. 5.2 as we found it to be necessary in order to tame the unbounded integrable singular surfaces that are of large dimensionality at these high loop counts 8 . The good agreement found for integral 3.6L4P.b.I is the first numerical confirmation of the analytical expression obtained in ref. [110].
Finally, the two entries 2L4P.a.I and 2L4P.b.I of table 2 present challenging integrals recently considered in ref. [111] (in which it appears as topology number B 72 ) in the context of the computation of the amplitude for Higgs production in association with a hard jet. In that work, the exact dependency on the internal quark mass is retained thanks to an original semi-numerical method for solving the system of differential equations relating master integrals. In the case of an internal top quark (2.2L4P.a.I), the authors could validate most of their results against the fully numerical ones obtained from sector decomposition techniques, however the case of the much lighter bottom quark (2.2L4P.b.I) proved to be more challenging for these approaches. The result from numerical LTD agrees with ref. [111] and has a numerical integration error only marginally impacted by the different values selected for the internal quark mass.
In fig. 20, we explore the stability of our numerical integration for two different classes of four-point kinematic configurations on one-, two-and three-loop ladder scalar integrals.
The first class of kinematics is unphysical, with p 2 1 = −5 and p 2 2 = p 2 3 = p 2 4 = (p 1 + p 2 ) 2 = −1. It is such that the region (p 1 + p 3 ) 2 = t > −7 can be addressed without any contour deformation, and for which we already showed results in fig. 1 of ref. [79]. In the complement region t < −7, a threshold singularity develops that corresponds to a single E-surface in this particular parametrisation of the kinematics and at any loop count. Our construction of the contour deformation involves 1, 4 and 11 sources (N S ) for the 1-  Fig. 20 demonstrated that numerical LTD is stable for different angular configurations, even when close to the crossing of thresholds in the external kinematics. We have however already observed in the one-loop results of table 1 that the convergence mostly depends on the shape and overlaps of the threshold singularity surfaces, which can become increasingly more complicated for boosted external momenta (that is | p i | 2 (p 0 i ) 2 ). In tables 4 to 8, we therefore seek to more systematically explore the performance of numerical LTD for external scattering 9 kinematic configurations p 1 p 2 → p 3 . . . p N of progressively stronger hierarchies in the scales m 2 j := p 2 j and s := (p 1 + p 2 ) 2 . We provide our explicit choice of kinematics in the ancillary material and we limit ourselves here to reporting their relevant scales: where the two different values for the masses of all internal propagators correspond to the massive (resp. massless) case labelled with (resp. without) a in the tables. We note that the series of kinematics K3 features internal propagators with masses set very slightly above that of one of the external momenta. This specific choice of internal mass is such that the existence condition of some E-surfaces are very close to being fulfilled, thus placing this challenging kinematic very close to crossing a threshold. Similarly to what can be observed in the scan shown in Fig. 20, we find numerical LTD to be in general stable even when approaching thresholds. At one loop (tables 4 and 5), we observe that the convergence mostly depends on the multiplicity of the external momenta, with a central value in agreement with MadLoop's reference beyond the percent level. At two loops (tables 6 and 7) and for integrals with more than four external legs, we find the scattering type of kinematics considered to be significantly more challenging than their decay counterpart featured in table 2 and we could not obtain a benchmark result from pySecDec. In those cases, the columns ∆[%] and ∆[%]| · | refer to the Monte-Carlo precision (and not the discrepancy w.r.t to the benchmark result) relative to the central value, and ∆[σ] is not applicable.
While numerical LTD generally performs well for kinematics featuring weaker hierarchies among its invariants, such as kinematics class K1, we found integrals where the convergence for the kinematics K2 and K3 was not good enough with our default deformation hyperparameters for the results to be reported in the tables. We note however that adjusting the two contour deformation hyperparameters th (which governs the strength of the expansion constraint), and M (which governs the strength of the anti-selection) can significantly improve the results. We illustrate this by optimising these two parameters for a particular six-point two-loop integral (2L6P.a) and for the K2 kinematics. Using a lowstatistics (50M points) exploratory scan, we find the optimal value of ( th ,M ) to be close to (0.7,0.01) for this configuration (most of the sensitivity lies in M ). We then report in the  table below  The two-loop eight-point integral 7.2L8P.K1 shows good convergence, but we could not obtain a result for its massive counterpart 7.2L8P.K1 * because it features a challenging maximal overlap structure (despite involving less than the 46 unique E-surfaces of integral 7.2L8P.K1) that we could not determine in a reasonable amount of computing time using the algorithm described in sect. 6.1. Beyond two loops (table 8), we again observe a significant improvement when considering massive internal propagators, which can partly be explained by the fact that in this case the deformation is no longer forced by the dynamic scaling of eq. (3.37) to become zero on the focal points of existing E-surfaces. We should mention that the four-point four-loop integrals included in the tables are at the upper end of the complexity that can currently be handled by our implementation. For massless internal propagators, the scattering kinematics Ki does not yield a good enough convergence while the decay kinematics necessitated an adjustment of the contour deformation hyperparameters (using a value for the parameter M in eq. 3.12 smaller than our defaults, e.g. M ∼ 0.01). Given that such integrals are also beyond what is of current phenomenological relevance, we present their results mostly to highlight the potential of numerical LTD.
Despite the wide range of variances obtained, we always find the central value obtained from numerical LTD integration to be within less than five sigmas away from the analytical benchmark ones (when available), as indicated by the ∆[σ] column of the tables. This observation is actually the most important aspect of our results, since in this work we first aim at demonstrating that our numerical implementation of LTD is robust and can therefore be predictive. Maximising numerical efficiency and exploring the optimisations discussed in sect. 5 is left to future work, for which results presented in this section can serve as a comparison baseline.

Divergent one-loop four-and five-point scalar integrals
We apply the subtraction scheme presented in sect. 4 to one-loop four-and five-point functions with massless propagators. For a randomly selected phase-space configuration, we go through all combinations of setting external momenta on-shell. For both the box and pentagon kinematics, we set s 12 = 1. For the box topology, when one of the external momenta is massive, we set m 2 1 = 1 4 , m 2 2 = 1 8 , m 2 3 = 2 9 , m 2 4 = 1 9 , respectively. For the pentagon topology, the masses are set to m 1 = 0.10, m 2 = 0.11, m 3 = 0.12, m 4 = 0.13, m 5 = 0.14. The results for these different configurations are shown in fig. 21, where the particular combination of masses for the external momenta is labelled by a binary number with the convention that a 1 in the ith position means that the ith external momentum is massless. We use the Cuhre integrator from Cuba package [100] with 200 million sample points. The time for each evaluation is independent of the mass configuration and is similar with the one presented in tab. 1.
Both the four-point ("box") and five-point ("pentagon") function can be integrated with high accuracy and precision: all but one of the central values are within a 0.005% of the analytical result. Only the imaginary part of the box topology with all the external momenta on-shell has a large uncertainty. The reason is that the central value of this integral is ten times smaller than for the other box configurations. However, even this point lies within 0.024% of the analytical result and has a relative standard error of 0.036%.
The analytic expression of the box integral and the triangle integrals required to construct the analytical expression for the counterterms have been computed using qcdloop [112]. The pentagon integral has been obtained using MadLoop5 [102,103] (ML5 henceforth).
7.3 One-loop amplitude for qq → γ 1 γ 2 and qq → γ 1 γ 2 γ 3 In this section we present the results from the integration of the amplitudes dd to two and three photons. For simplicity, we kept the order of the final photons fixed; the actual result for the amplitude can then be recovered by permuting through the final-state photon momenta. The helicities are defined following the HELAS convention [113], and are taken positive for all the external particles. The evaluation of the numerator, involving contractions of Lorentz and spinor indices, is performed numerically at run-time. This is not an efficient way to perform the numerator algebra, but the aim of this work is to highlight how LTD can be used to obtain results for physical and divergent expressions.
The analytic expressions have been compared with ML5 with g s = 1.21771, g = 0.30795 and µ r = 91.1880 as couplings. We also remind the reader that the results from ML5 are rescaled by an overall factor (4π) /Γ(1 − ).
For the dd → γ 1 γ 2 process, we consider the process in its centre-of-mass rest frame, with the quarks aligned along the z-axis. The result will only depend on the scattering energy and angle. The former is kept fixed and corresponds to a simple rescaling of the integral and the latter is varied in a scan and plotted in fig. 22. We used the Cuhre integrator from Cuba package [100] with two million evaluations. In the last plot of fig. 22 we notice that the result is almost completely determined by the integrated counterterms. This is especially true for the real part, where one can see that resulting regulated integral is six orders of magnitude smaller than the finite part of the analytic expression.
As for the case of scalar divergent integrals, we use the Cuhre integrator with however only 2 million sample points in this case. Despite this relatively low statistics, a large fraction of the results already have relative error below 0.05%. In the upper plot of fig. 22 we show the relative deviation with a large scale in order to highlight the few points that are not within this small error. One important observation however is that the Monte-Carlo error reported is reliable, as highlighted by the fact that all discrepancies are smaller than one (in modulus) when expressed in unit of the Monte-Calo standard deviation σ.
In fig. 23 we show a scan of dd → γ 1 γ 2 γ 3 . In the same way as for the two-photon production case, we consider the scattering in the centre-of-mass rest frame. This time however, the number of unspecified and non-trivial degrees of freedom is four so that keeping a fixed energy s 12 = 1 leaves us with three parameters. For the kinematic configuration d(p 1 )d(p 2 ) → γ 1 (p 3 )γ 2 (p 4 )γ 3 (p 5 ), we choose to scan in the angle θ 13 = ∠(p 1 , p 3 ), and s 45 which gives an indication of how collinear the momenta p 4 and p 5 are. We fix the remaining degree of freedom by forcing the process on a plane, which allows for the configuration where p 4 is collinear to p 1 , thus resulting in the valley shown in plots (a -b) of fig. 23. For dd → γ 1 γ 2 γ 3 , we observe that the relative contribution from the integrated counterterms is not as large as for dd → γ 1 γ 2 , because this five-point amplitude has more contributions that are IR-finite (specifically D4-6 from fig. 15) and therefore not captured by the counterterms.
We can see that the relative error is < 1% for most of the points in the scan as shown in the upper plane of plot (e-f) from fig. 15). In the lower part of the same plots, the precision of the result with an error that is also < 1% for most of the points. Along the valley, the relative accuracy is not as good as in the other regions, which is to be expected when the central value of the integrated expression becomes smaller that the values around it. Similarly as to elsewhere in this subsection, the results were obtained using the Cuhre integrator and 2 million sample points. The low number of samples is due to two mainly two reasons: first, we used a naive implementation of the numerator containing spinor chains that are recomputed numerically for each evaluation and second, despite the measure taken for improving the UV behaviour of the integrand, probing that region still requires many evaluations in quadruple precision thus increasing the overall evaluation time by roughly one order of magnitude when compared to the corresponding scalar topologies. In the present work, we put no effort in optimising the numerator expression which we leave to future work. The main objective of these results is to demonstrate the viability of computing physical amplitudes with numerical LTD by combining the contour deformation together with the necessary infrared and ultraviolet counterterms. Optimising the implementation of the numerator will allow us to handle more complicated amplitudes and to consider higher integration statistics.    [112]. The (nominal) horizontal axis shows different phase-space configurations using a binary notation, where a 1 (resp. 0) in the ith position signifies that the ith external momentum is on-shell with p 2 i = 0 (resp. off-shell, that is with p 2 i = 0). All but one of the central values are within 0.005% of the analytical result. The outlier with configuration 1111 lies within 0.024% of the analytical result and has a relative standard error of 0.036%.    -73 -

Conclusion
The ongoing and future research programme of the LHC calls for improving on the theoretical accuracy of the simulation of many scattering processes. A formidable effort from the high energy physics community over the last decades lead to the computation of many higher-order corrections of key relevance. However, computing QCD amplitudes beyond two loops and/or four scales remains extremely challenging, even with modern analytical techniques. We identify this problem as being one of the main bottlenecks whose resolution demands a radically new approach.This observation is what motivates our work on numerical Loop-Tree Duality, as its strength and limitations are orthogonal, and thus complementary, to those of the canonical paradigms for predicting collider observables. The potential of numerical LTD is reinforced by the promising perspective it entails regarding its eventual combination with real-emission contributions. In our recent work of ref. [79], we presented our first developments and generalisation of LTD and, encouraged by our findings, we proceeded in this work to extend its range of applicability. First, we established a contour deformation for regulating the threshold singularities exhibited by loop integrals when considering physical scattering kinematics. In accordance with our long-term goals, we built a solution that is prone to automation and made no compromise regarding the generality of numerical LTD: availability of computational resources should remain the only limiting factor. Moreover, we insisted that the validity of the contour deformation should be independent of the particular values of its hyperparameters, thus guaranteeing the predictive power of numerical LTD. We demonstrated that our construction and implementation achieves these objectives by applying it to over 100 different representative configurations, ranging from one-loop boxes to four-loop 2x2 fishnets.
Second, we presented our first step towards computing divergent integrals and physical amplitudes. This requires combining the LTD expression with local integrand-level counterterms regularising divergences occurring for ultraviolet, soft and/or collinear loop momenta configurations. We described this subtraction procedure at one loop and showcase explicit examples for divergent scalar four-and five-point integrals, as well as for the oneloop amplitude of the production of two and three photons. This paves the way for a first application of numerical LTD to the numerical computation of two-loop divergent scalar integrals and of complete two-loop amplitudes, using the local counterterms introduced in refs. [81,82].
In this work, we focused on further developing numerical LTD in a way that is provably correct, general and that demonstrates predictive power. Therefore, we did not tune our hyperparameters for the hundreds of cases we studied and, although already satisfactory, the numerical convergence and run-time speed showcased by our results are by no means final. We leave their improvement to future work.
The ability to locally regulate ultraviolet and infrared singularities at higher loops and the performance of the numerical convergence are two key difficulties whose resolution will determine the eventual viability of numerical LTD. Our work shows a clear path for this novel approach to significantly contribute to the effort of meeting the theoretical accuracy goal set by the needs of current collider experiments.
-74 -We would like to thank Francesco Moriello for providing us with the opportunity of applying our method to a case of practical interest, Stephen Jones for helping us produce comparison results with pySecDec, Mao Zeng and Babis Anastasiou for providing insights on the amplitude subtraction and Armin Schweitzer and Rayan Haindl for fruitful discussions. This project has received funding from the European Research Council (ERC) under grant agreement No 694712 (PertQCD) and SNSF grant No 179016. Numerical results presented in this work used computational resources from the Piz Daint cluster, administered by the Swiss National Supercomputing Centre (CSCS).