Direct numerical evaluation of multi-loop integrals without contour deformation

We propose a method for computing numerically integrals defined via iϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \epsilon $$\end{document} deformations acting on single-pole singularities. We achieve this without an explicit analytic contour deformation. Our solution is then used to produce precise Monte Carlo estimates of multi-scale multi-loop integrals directly in Minkowski space. We corroborate the validity of our strategy by presenting several examples ranging from one to three loops. When used in connection with four-dimensional regularization techniques, our treatment can be extended to ultraviolet and infrared divergent integrals.


Introduction
The ever-increasing precision of data from particle physics experiments requires a comparable or better level of precision in theoretical predictions, both to establish the parameters of the Standard Model and to search for physics beyond it. To achieve such precision requires the computation of multiloop amplitudes. A fundamental ingredient of such calculations is the evaluation of master loop integrals (MIs), in terms of which the problem is reduced. This can be performed by analytic, semi-numerical or fully numerical techniques (see [1] for a recent review). Analytic methods are very successful when the class of functions that contribute to the result is known, which usually happens when the number of internal and external masses is limited. However, such a-priori knowledge is not always available, especially when the number of scales increases, so that in these cases one would like to be able to compute MIs numerically, for instance by Monte Carlo (MC) techniques.
In the numerical computation of MIs, an important problem is the appearance of integrable threshold singularities, a e-mail: pittau@ugr.es (corresponding author) b e-mail: webber@hep.phy.cam.ac.uk where single poles are moved away from the real integration domain by the i prescription. These singularities require special treatment, such as a contour deformation into the complex plane [2][3][4][5], or vanishing-width extrapolations methods [6][7][8][9]. Contour deformations are usually controlled by some parameter whose value should be not too small, to guarantee numerical accuracy, and not too large, to avoid crossing branch cuts. In extrapolation methods a series of integrals should be determined that converges to the right value while keeping the computation time low.
This paper explains how integrals defined through the i prescription acting on first-order poles can be evaluated numerically without deforming the integration contour into the complex plane, and how this can be employed to compute MIs appearing in multi-loop calculations. In addition, we demonstrate that this strategy allows one to compute recursively higher-loop functions in terms of lower-loop ones. Other semi-numerical methods relying on one-looplike objects to build higher loops can be found in [10][11][12]. In [10] a Wick rotation of the loop momentum is needed to avoid singularities. The Feynman parameter space is used in [11], and a contour deformation in [12]. Our method works directly in Minkowski space and avoids contour deformations.
The structure of the paper is as follow. Section 2 details our approach. In Sect. 3 we use it to integrate numerically threshold singularities after an analytic integration over the energycomponents of the loop momenta. Section 4 explains how to glue together lower-loop structures to compute numerically certain classes of higher-loop MIs. Finally, in Sects. 5 and 6 we extend our treatment to ultraviolet and infrared divergent configurations regularized via the four-dimensional method of [13].

Avoiding contour deformation
In this section we present two methods which avoid contour deformation. The first method uses complex analysis, while the second approach directly works with the original integrand. The two procedures are equivalent, in that they give rise to the same mappings. The numerical results presented in the paper are obtained with method 1 and cross-checked with method 2.

Method 1
For the sake of clarity, distinct letters (with or without additional subscripts) are used to denote variables ranging in different intervals. In particular, we employ x when −1 ≤ x ≤ 1, y if 0 ≤ y ≤ 1, σ provided −∞ < σ < ∞. Finally, 0 ≤ ρ ≤ 1 stands for a random Monte Carlo (MC) variable.
The core of the procedure is a change of variable such that the 1/(x + i ) behaviour of the integral is flattened with x ∈ R. This is obtained by imposing where z is a new complex integration variable. In fact, inserting (2) in (1) gives the desired result, Equation (3) evaluates to −iπ along any curve in the z complex plane connecting z = 0 to z = 1 when → 0. We use this freedom to impose x ∈ R by parametrizing z = α + iβ with α, β ∈ R and π ≤ α ≤ 1 − π .
Inserting (4) in (2) gives which is real when πβ = ln sin[π(1 − α)] , namely Therefore which gives where g has been introduced to impose the normalization to 1 also for small but not vanishing values of . In summary, after changing variable as in (2), the requirement x ∈ R determines the relation between e(z) and m(z). Armed with these results, we generalize (1) to an integration over a function with φ(x) sufficiently smooth at x = 0, 1 Splitting the integration region of (8) into the two sectors with x α < 0 or x α > 0 gives (11) where y α := /tan(απ ). Equation (11) can be translated to a MC language by looking for the local density g(y) that corresponds to a change of variable dρ = g(y)dy reabsorbing the singular behaviour of the integrand of (10), with´1 0 dy g(y) = 1. By comparing (12) to (11) one determines The mapping of (13) optimizes the integration over the real part of z (see (7)). This gives stable numerical results when φ(x) is such that the y α / terms in (11) are suppressed. When this is not the case, they generate a large contribution to the variance and, in order to flatten them, the parametrization complementary to (7) is necessary, which gives where Again, (15) is correctly normalized also for a small but not vanishing . Comparing (12) to (15) gives now
To reduce the variance when φ(x) peaks inside −1 ≤ x ≤ 1 in a known way, it is also possible to include an arbitrary numbers of further channels g i (y) (i > 2). However, care must be taken due to the fact that the MC weight of (12) includes both the f (−y) and f (y) contributions. To determine the corresponding density we observe that which means that if x is randomly chosen in −1 ≤ x ≤ 1, the density is g i (−|x|) + g i (|x|). Hence, the MC weight is f (−y)+ f (y) / g i (−y)+ g i (y) , with g i normalized such that In summary, with N ch channels (including g 1 and g 2 ), the more general multichannel MC mapping reads where with arbitrary (but self-adjustable) weights fulfilling In the actual MC used to produce the results presented in this paper we superimpose on g c (y) a flat distribution g 3 (x) = 1/2 and a channel g 4 (x) = 1 2 ln 1+δ which takes care of peaks around |x| = 1.

Principal value integrals
It is often useful to deal with improper integrals, whose behaviour at large values of the integration variables is defined via the Cauchy principal value. The fact that the two symmetric points with respect to x = 0 are always considered together, makes the use of (20) very convenient. As a matter of notation, we define which can be mapped onto the interval [−1, 1] by changing variable, Thus, for instance, where we understand the symmetric treatment of (20), so that (25) is well defined even when φ(σ ) approaches a constant as σ → ±∞.

Multiple integrals
Equation (20) can be easily extended to n-fold integrals of the type Our notation is such that {x} = x 1 , x 2 , . . . , x n and φ({x}) is a smooth function at {x} = {0}. The result is where dρ j = g tot (y j )dy j and the numerator stands for a sum over the 2 n terms with positive or negative arguments. For instance, when {y} = y 1 , y 2 , Equation (27) can be generalized to more poles per variable, moved away from arbitrary domains ∈ R, either by partial fractioning the integrand or by splitting the integration region into sub-intervals. However, configurations like that never appear in what follows, so we do not pursue a detailed analysis in this direction.

Method 2
As an alternative to the above method, one can apply separate changes of variables to flatten the real and imaginary parts of the pole factor(s) in the integrand. Consider the integral . We can write this as where r m := ln(1 + a 2 / 2 )/2 , θ m := arctan(a/ ).
Each of these integrals has optimal variance reduction (in the absence of information about φ) and is therefore suited to numerical integration as long as φ(x) is smooth at x = 0. Note that the /(x 2 + 2 ) and x/(x 2 + 2 ) behaviours of (29) correspond to the local densities in (13) and (16), respectively. The method is easily generalised to two variables. Consider where For MC evaluation, we proceed as follows: for each shot, where 0 < r 1,2 < r m , 0 < θ 1,2 < θ m uniformly, with r m and θ m as in (31). In (34), set x 1 = x 1t when the first superscript is 0 and x 1 = x 1r when it is 1, and similarly for x 2 according to the second superscript. The weights for the real and imaginary parts are then The generalisation of (32) to n variables is clear: φ {k j } has superscript k j = 1 in the jth location when there is an x j in the integrand, otherwise k j = 0. The symmetrized function becomes where For MC evaluation, for each shot, generate two points in the n-dimensional hypercube where again 0 < r j < r m and 0 < θ j < θ m uniformly. In (39), set x j = x jt when k j = 0 and x j = x jr when k j = 1. The weight is then Note that each shot involves 2n random numbers for {x jt } and {x jr }, and then 4 n function evaluations at x j = ±x jt and ±x jr , so the computation time increases rapidly with the number of variables.

Choosing
Here we perform a study of the value of to be used in practice. More specifically, we compare the numerical and analytic determinations of the three-fold test integral whose behaviour at x j ∼ 0 mimics a typical multidimensional environment. The result of this comparison is given in Fig. 1 Requiring the = 0 bias on Q MC ( ) to be of the order of the maximum between and the relative MC error gives the condition Table 1 reports the R value of the entries of Fig. 1 and their MC accuracy defined as From Fig. 1 and Table 1 we infer that a range 10 −8 ≤ ≤ 10 −6 is adequate to achieve MC estimates accurate at the level of three parts in 10 5 . Since the results presented in this paper are never more accurate than this, we set, for definiteness, = 10 −7 . However, the last row of Table 1 shows that numerically stable predictions are produced also with a smaller and a larger MC statistics. From this, we deduce that the = 0 bias can be reduced to be negligible in most practical applications, and that the accuracy of our method is driven by the MC error.

A semi-numerical integration algorithm for MIs
In this section we illustrate how the approach of Sect. 2 can be successfully applied to produce stable and precise seminumerical MC estimates of loop MIs. This is achieved in  two steps. Firstly, we integrate analytically over the energy components of the loop momenta, which is always doable by means of the Cauchy integral theorem. In addition, depending on the case at hand, some of the loop angular integrals can also be performed analytically. In this way, integral representations of MIs can be easily obtained. Secondly, we give up any attempt towards a fully analytic integration, which may be difficult, and integrate numerically over the leftover loop components. The integrand to be evaluated is usually plagued by threshold singularities. Single poles migrate towards the real integration domain for some kinematic configurations, so that a blind numerical integration over denominators deformed by the Feynman i prescription gives large errors. However, this is precisely the situation for which our approach is designed. We mitigate these problems by retaining a finite small value of , flattening the real and imaginary parts of pole contributions, and applying multichannel mappings. 2 In what follows we illustrate the performance of this strategy by means of two examples. gives with One splits where R ± 2 := ( The cut of the logarithms with +i (−i ) is in the lower (upper) t complex half-plane, so that the integration over t in (49) is trivial once one rewrites The results is where When √ τ > 2, the first integrand of (52) develops a pole at r = i that migrates towards the integration region in the limit → 0. Treating this with the strategy of Sect. 2 gives the results presented in Table 2.

A two-loop example
We study the two-loop self-energy scalar diagram of Fig. 3. For a timelike P/m = ( √ τ , 0) and m 1 = m it reads where and with μ 0 := m 2 0 /m 2 . Integrating over the angular variables gives The integration over t 1 and t 2 is trivial and produces where When m 0 = 0, a two-dimensional implementation of the method of Sect. 2 gives the results reported in Table 3. Larger MC errors correspond to smaller values of ρ. However, we observe that when μ 0 = 0 this effect is mitigated. For instance, a 10 9 MC-point estimate with ρ = .1 gives

Gluing together lower-loop structures
Here we show how higher-loop integrals can be expressed in terms of lower-loop building blocks. Throughout this section dimensionful quantities are rescaled by an arbitrary mass m, so that loop momenta are written as in (54) and, in particular, Furthermore, we define Table 3 The two-loop integral (57) with m 0 = 0 multiplied by −m 2 τ/π 4 for several values of ρ := 4/τ . Numbers obtained with 10 9 (10 10 ) MC points when ρ > 1 (ρ < 1). The analytic result is taken from [17]. MC and study cases up to a P := p 1 + p 2 → p 3 + p 4 kinematics of the form Rescaled propagators belonging to the loop momentum q are denoted by In addition, we define The essence of the procedure is to use σ 0 and σ 1 as integration variables of the method of Sect. 2. This is achieved by multiplying the integrand by where This gives rise to the appearance of the following three functionals, where j = 2, 3. Assuming J 1 independent of any angular variable and J 2 (J 3 ) independent of θ (φ) allows one to compute the functionals once for all. 1 reads As for 2 , one has with Note that (66) and (67) have been analytically continued to configurations with any sign of τ , τ 1 , τ 2 , λ 12 . Finally in which  67) and (69) produce In (71) the integration over the azimuth angle φ is trivial, while the integral over c θ in (72) has to be dealt with numerically using the method of Sect. 2, which gives where

Improving the numerical accuracy
When inserted in (71) and (72) Note that case (a) is relevant to C but not to D 0 , since while (b) applies to both C and D 0 due to the common 1/τ prefactor. In the following paragraphs we illustrate how numerical inaccuracies caused by the configurations (a) and (b) in (76) can be circumvented.
As for case (a), a preliminary analysis is in order to understand the mechanism that makes C finite despite (77). 4 We define in terms of which Now the σ 0 integral is convergent by power counting and λ in (61) behaves as where α and β are constants. Replacing λ withλ in (73) produces an integrand in which all branch points and poles are located in the lower σ 0 complex half-plane. As a result, the integral over σ 0 approaches zero when σ 10 → ±∞, so that the → ∞ limit exists. This same reasoning allows one to construct a class of vanishing integrals defined as with m 2σ 4 The principal value integrations in (71) are not sufficient to regularize the large σ 0,1 behaviour. In fact, (μ0,μ1) 2 [1] approaches two different constants when σ 0,1 → ∞ or σ 0,1 → −∞, so that no cancellation is possible.
whereq is the asymptotic |σ 10 | → ∞ limit of q, Again, all cuts and poles lie in the lower σ 0 complex halfplane, so that Now α and β can be set to obtain a local cancellation of the problematic large σ 10 configurations. 5 An explicit calculation with τ 1 = τ 2 = 0 gives When τ i = 0 the accuracy of (71) is improved by the nonvanishing external masses, so that (85) is relevant to this case as well. In summary, the formula produces numerically stable results with α and β given in (85) when the same sequence of σ 0 and σ 1 values are used in both C andC. It turns out that the configurations of type (b) of C are also cured by the subtraction in (86). Thus, we are only left with the discussion of the case (b) for D 0 . In the τ → 0 region it is convenient to give up the exact formula and use, instead, a few terms of a Taylor expansion in τ and χ obtained with the method given in Appendix A, By doing that, it is easy to find a value of τ below which the exact result is well approximated by (87), and above which (72) is accurate.
The results for C(P 2 , 0, 0, m, m, m) are shown in Figs. 5 and 6. In the latter, the relative difference between the MC and the analytic (AN) non-zero results of the former is plotted in terms of where R and I refer to the real and imaginary parts, respectively. With the given MC statistics, the analytic result is  [16]. Results obtained with 10 9 MC shots per point reproduced by (86) within a few parts in 10 4 , 6 which is an accuracy comparable to the one of Table 2, but obtained with 50 times more points. The main reason for this difference is that the one-dimensional representation (52) is now replaced by the twofold integration (71). On the other hand, the gluing algorithm is highly modular and can be extended to more complex situations, as we will see in the next subsection. Table 4 displays our results for D 0 as a function of τ and the scattering angle −χ/τ = (1 + cos θ 13 ) /2. We use the MC evaluation of (72) above τ = 1 and the Taylor expansion of (87) when τ < 1. In the former case, an accuracy of a few parts in 10 3 is reached with 10 9 MC points.

Two-loop self-energy
The two-loop diagram of Fig. 3 can be easily obtained by gluing together a one-loop triangle and a tree-level structure. For instance, with m 0 = m 1 = 0 one has Equation (89) suffers the inaccuracies of type (a) of (76). To cure this, we use the same strategy described in Sect. 4.1. First, we consider an explicit representation of the triangle as a function of σ 0 and σ 10 in (79), with λ 0 in (75), and From this, it is easy to determine a |σ 10 | → ∞ asymptotic approximant of (90) that gives zero upon integration over σ 0,1 . The results reads 7 whereλ 0 := (σ 10 −τ ) 2 −4τ , c is arbitrary andF(σ 0 , σ 10 ) is constructed by replacing in (91) the u j , v j with their asymptotic counterpartsũ j ,ṽ j defined as Finally, in the same spirit as (86), we rewrite where we understand a local subtraction of the large |σ 10 | configurations.
In Table 5 we present our numerical estimates based on (95) with c = 1/2. An accuracy of the order of 10 −3 is achieved with 10 9 MC points. We also studied the stability of (95) at small values of τ . For instance, when ρ = 100 (τ = 0.04) we obtain Note that determiningS 2 ( c ) requires an analytic knowledge of the integrand. When this is not possible, an alternative approach is to cut away the problematic configurations in a controlled manner. For instance, discarding in (89) integration points with is expected to produce an error of O( 2 ). Indeed we checked that, with = 0.01, the fraction of the integral discarded by (97) is always below the errors reported in Table 5.

Three-loop self-energy
Our next example is the scalar three-loop self-energy of Fig. 7, which is attained by gluing together the two trian- Table 5 The two-loop integral (95)  gles C L := C(τ, σ 1 + μ 1 , σ 0 + μ 0 , μ 2 , μ 3 , μ 4 ), by means of 1 , Now the integrand vanishes fast enough at large σ 0,1 , so that no subtraction is needed. Our results for the case {μ k } = {1, 1, 0, 0, 0, 0, 0, 0} are shown in Fig. 8, where we compare them with digitized curves obtained from Fig. 4 of [10]. The agreement is good, but our points tend to overshoot the lines at large τ . However, the quality of the plot in Ref. [10] is poor there and the digitization may be misleading. Thus, we cross-checked internally our high-energy results by comparing the outcome of two independent MCs based on method 1 and 2 of Sect. 2, respectively. We did

Planar two-loop vertex
Consider now the two-loop scalar vertex of Fig. 9. Gluing the one-loop triangle on the left by means of 2 gives This representation holds true with any sign of τ, τ 1 , τ 2 and for any choice of internal masses. 8 In addition, it does not require subtracting large σ 0,1 configurations. In Table 6 we collect a few results obtained by computing the triangle with OneLOop. The last row refers to the Standard-Model-like Z → νν case with m 0 = m 1 = m e , m 2 = m 3 = m 4 = M W , m 5 = 0, P 2 = M 2 Z and p 2 1,2 = 0. This shows that the MC error is under control also for configurations with large mass gaps. Note that V 2 ({μ k }) can also be obtained by gluing together the box on the right D R := D(τ 1 , τ 2 , σ 3 + μ 3 , σ 4 + μ 4 , τ, σ 34 b , μ 1 , μ 2 , μ 0 , μ 5 ) and the tree-level decay on the left. When k 2 is greater than zero one has The loop momentum of the m 3 line of Fig. 9 flows through just one propagator of D R . Because of that, a "technical" cut τ 2 /(τ 2 + (σ 3 − σ 4 ) 2 ) < is required to damp the inaccurate large σ 3,4 behaviour of (101). By power counting, the discarded contribution is of O( 4 ). With = 0.05 and 10 9 MC points (101) reproduces the numbers of Table 6 at the percent level.

Non-planar two-loop vertex
The same reasoning leading to (101) allows one to compute the non-planar two-loop vertex of Fig. 10, where D np := D(τ 1 , σ 3 + μ 3 , τ 2 , σ 4 +μ 4 , σ 34 a , σ 34 b , μ 1 , μ 2 , μ 0 , μ 5 ).  Now the loop momentum of the m 3 line flows through two propagators of D np . This provides an additional damping factor in (102) with respect to (101), so that large σ 3,4 configurations do not lead to numerical inaccuracies and no technical cut is required. A few numerical results are collected in Table 7.

Planar and non-planar double box
The planar and non-planar two-loop double boxes are depicted in Fig. 11a, b, respectively. They read  Fig. 11a, b, respectively. In Figs. 12 and 13 we present a comparison between our estimates and the results presented in [7] for the case in the region of τ where | cos θ 13 | ≤ 1. 9 The agreement is very good. As benchmark values, we list in Table 8 the MC entries of Figs. 12 and 13, together with their statistical errors.  [7]). Red bullets (blue squares) refer to the real (imaginary) part obtained with 10 9 MC shots per point using OneLOop to compute D R . The solid and dashed lines are derived from Fig. 3 of [7] In all cases presented so far we could perform an easy analytic integration over at least one angular variable of the loop momentum ω. When the number of loops and legs increases, integrating analytically over c θ and/or φ is not trivial any more, so that the number of numerical integrations required by the gluing procedure reaches its maximum value, i.e. four. In this paragraph we give a few examples of the gluing approach in these more complex situations. In particular, we use the multichannel approach of Sect. 2.1 to study the scalar triple box B 3 and two-loop pentabox E 2 of Figs. 14 and 15, respectively.
Note that it is not difficult to deal with non-planar configurations. For instance, the diagram obtained from Fig. 14 by interchanging the vertices v 1 ↔ v 2 and v 3 ↔ v 4 is as in (106), but with Finally, the pentabox depicted in Fig. 15 reads where E(ω) is the one-loop pentagon on the right side, which depends upon ten independent invariants -built from the momenta q = mω, p 3÷5 -and the five masses m 3÷7 . The presence of the 1/(σ 2 +i ) propagator forces one to trade the integral over c θ for an integration over σ 2 to be dealt with the method of Sect. 2. This, together with the change of variable of (108), gives where To provide a benchmark value, we have chosen to evaluate (113) with μ 0÷2 = 1, μ 3÷7 = 0 at a particular p 1 + p 2 → p 3 + p 4 + p 5 phase-space point satisfying τ = 50, τ 1 = τ 2 = 0 and τ 3÷5 = 1, namely 1, 0, 0) , Fig. 15 The two-loop planar pentabox E 2 of (112) r, r, 0) , with r = √ 193/450, c 34 = −159/193 and s 34 = 1 − c 2 34 . We employed CutTools [18] to reduce E(ω n ) to one-loop boxes. Computing the latter with OneLOop gives, with 10 9 MC shots, Once again, non-planar configurations are obtained without extra effort. For instance, if the vertex v of Fig. 15 is moved to the m 3 line, (112) still holds by simply modifying the oneloop pentagon E(ω) accordingly.

UV divergences
We deal with ultraviolet divergent integrals following the FDR approach of [13], in which the divergent configurations are extracted from the original integrand by partial fractioning. 10 The resulting expressions are integrable in four dimensions and nicely match the algorithm of Sect. 2. We illustrate our procedure by means of the scalar two-point function of Fig. 16, wherē with M 2 1 (q) := m 2 1 + 2(q · P) − P 2 andq 2 := q 2 − μ 2 + i . By partial fractioning 11 This gives, by definition of FDR integration, where μ R is the finite renormalization scale. It is convenient to take the limit μ → 0 directly at the integrand level and substitute μ with μ R only in the logarithms. This is achieved by rewriting where D i :=D i + μ 2 . By doing so, μ 2 is dropped everywhere, except in the logarithmically UV divergent part of the vacuum, where it is replaced by μ 2 R . In this way, no μ → 0 limit is required, so that (120) is a good starting point for a numerical treatment. In what follows we describe how the methods of Sects. 3 and 4 can be adapted to deal with (120). More complex multi-loop UV divergent configurations can be treated likewise.

Integrating over the loop energy component
We take m 1 = m 0 = m for simplicity. A rescaling as in (47) produces where with The Cauchy integral theorem allows one to compute Inserting this in (121) and using r = R 0 − √ τ /2 as a new integration variable gives where If √ τ > 2 the pole at r = i migrates towards the integration contour when → 0. Treating this with our numerical approach produces the results collected in Table 9.
Integrals with a polynomial degree of divergence can be treated in exactly the same way. As an example, Appendix B details the case of the one-point function

Gluing substructures
The gluing approach of Sect. 4 can be easily extended to (116). Inserting (63) in (120) gives with and ν defined in (123). The presence of a double pole is an obstacle to a direct numerical treatment of (128). In fact, our algorithm is designed to deal with single poles only. However, we observe that, if ν has a finite imaginary part, the singularity never approaches the real axis. In particular, B FDR (−iν) is better suited than B FDR (ν) to be evaluated numerically if ν ∈ . Besides, the connection between the two can be derived by differentiating (120), which gives B FDR (−iν ) still suffers from numerical inaccuracies of type (76) (a). To cure this, we locally subtract from it an approximant,B FDR (−iν ) = 0, constructed in such a way that, after changing variables as in (79), all cuts and poles lie in the lower σ 0 complex half-plane. This is obtained by replacing  Table 10 The integral of (131) with ν = 1/4 and μ 0 = μ 1 = 1. The estimates are obtained by sampling with 10 8 MC shots (132) evaluated at ν = 1. MC errors between parentheses In Table 10 we present our estimates for B FDR (1/4) with μ 0 = μ 1 = 1 obtained by means of Eqs. (131) and (132). The figures match the results of Table 9, although with larger errors. However, we point out that the gluing method is more flexible when it comes to generic kinematics. For instance, with τ = −10, μ 0 = 1, μ 1 = 4, one obtains, with 10 8 MC points, to be compared to the analytic value −i 27.220.

IR divergences
We deal with IR divergent integrals by means of the FDR approach of [19], where a small mass μ, added to judiciously chosen propagators, is used as a regulator of both infrared and collinear divergences. In this section, we illustrate how this allows one to combine virtual and real contributions prior to integration. After that, our method can be used to evaluate numerically loop integrals where the IR configurations are locally subtracted. We study, in particular, the IR divergent scalar triangle that appears in a P → p 1 + p 2 decay with p 2 i = 0 and s := P 2 . However, our findings can be generalized to more complex environments.
Our strategy is based on combining together cut-diagrams that are individually divergent, but whose sum is finite. We use a scalar massless gϕ 3 theory defined through the Feynman rules of Fig. 17, where we have introduced propagators with positive and negative values of the energy and the momentum component along the x direction. The cuts contributing to ϕ * → ϕϕ(ϕ) are listed in Fig. 18 where, to make contact with (134), D a + D e = 4π 2 g 4 π 2 (i C IR ), The diagrams are organized in pairs sharing collinear singularities. For instance, in D a the energy component of propagator 1 is the sum of those of propagators 2 and 3. Thus, propagators 1 and 3 never pinch in the q 0 1 complex plane, and propagator 3 can only become collinear to 4. Likewise, in D b Fig. 17 The Feynman rules of the gϕ 3 theory. A special notation is used for propagators with positive and negative values of p 0 and p x , which, in our convention, coincides with the direction of the back-toback final-state particles in the P rest-frame. The complex conjugate of such rules is used in the r.h.s. of diagrams cut by a dashed line Fig. 18 The two-and three-particle cuts contributing to ϕ * → ϕϕ(ϕ) in the P rest frame where the propagators 2 and 5 have negative and positive components of the momentum along x, as indicated by the ∓ labels attached to them. Thick lines represent the propagators to which μ 2 → 0 is added, as explained in the text the sign of the momentum components along x only allows particles 3 and 4 to become collinear to 5. In both cases, we regulate the singular splitting by including a small mass μ in propagators 3 and 4, leaving 5 massless. 12 In summary, D a + D b is free of collinear divergences, and the same happens for D c + D d . In addition, D a + D b + D c + D d is also free of infrared singularities. A similar reasoning applies to D e + D f + D g + D h , but with an opposite sign of the energy component of q 1 . The previous analysis shows that the threeparticle cuts D b and D d can be used as local countertems for D a + D c . 13 This requires common reference frames. One can employ two different routings for D a +D b and D c +D c . However, they must coincide in the limit q 1 → 0 to guarantee the cancellation of the soft behaviour of D a + D b + D c + D d . In particular, when computing D a,b we assign a momentum q 2 to propagator 4, from left to right, and choose On the other hand, we calculate D c,d with q 2 assigned to propagator 5 and The result of the computation is reported in Appendix C in terms of integrals over It is convenient to further split D a,c = D s a,c + D u a,c , where the superscripts s,u refer to the subtracted and unsubtracted regions, which correspond to the integration intervals √ η < R 1 < 1/2 and 1/2 < R 1 < ∞, respectively. In fact, D b,d contribute in the subtracted region only, and D u a + D u c is free of IR singularities, s g 4 (D u a + D u c ) = −2π 5 ln 2 2 + 2Li 2 (−1/2) = 254.838137 · · · .
An analytic calculation [19] shows that j=a,b,c,d D j = 0. Hence, one must have In Table 11 we display our numerical estimate of K based on Eqs. (C.22), (C.29) and (C.34). The correct result is precisely approached and the MC error does not grow when decreasing η, which is an indication that the local cancellation works as expected. Finally, we point out that the outlined strategy can be turned into a fully exclusive local subtraction algorithm by introducing suitable phase-space mappings, as described in [20].

Conclusion and outlook
We have presented a flexible method for the numerical treatment of loop integrals in four-dimensional Minkowski space, without the need of explicit contour deformation. This is achieved by exploiting the i prescription with a small finite value of and making changes of variables to reduce the variance of both the real and imaginary parts of the integrand. We propose a semi-numerical approach, in which an analytic integration over loop time-components is followed by multichannel Monte Carlo integration. In some cases, further integrations can be performed before the final numerical step. The method lends itself readily to the evaluation of complex multi-loop structures by gluing together simpler substructures. It also deals easily with processes involving many different external and propagator mass scales, where analytical results are difficult to obtain. In practice, we find that 10 9 Monte Carlo shots with ∼ 10 −7 (in terms of some relevant mass scale) can yield relative precision of the order of 10 −4 for one-loop diagrams and 10 −3 for two-and three-loops obtained by gluing together analytical results for one-loop substructures. As for the performance of our algorithms, we report in Table 12 the time to produce 10 6 MC shots with method 1 for a few representative cases. It ranges from a few tenths of a second to more than a minute. Method 2 gives somewhat slower timings.
We have focused on scalar integrals without any structure in the numerator, but we expect that the treatment of loop tensors should follow the same guidelines described in this paper. In particular, the approach of Sect. 3, in which an analytic integration is performed over the loop time-component, should work as it stands. As for the gluing method of Sect. 4, adding structures in the numerator could potentially lead to worse behaviour that needs to be corrected by local subtractions of large loop configurations, as done in Eqs. (86), (95), (132), or by the technical cuts described in (97) and (101). We leave a detailed study of this subject for further investigation.
We have sketched out how our method can be extended to UV and IR divergent configurations. Again, a deeper investigation is left for the future.
In summary, we believe that a numerical treatment of virtual corrections in four dimensions, of the type we have proposed, could be very beneficial in the computation of complicated multi-leg multi-scale amplitudes. More specifically, we think that the direction to go would be to integrate directly the amplitude as a whole, rather than the separate MIs. This could mitigate some of the large gauge cancellations among individual contributions, if common loop momentum routings are chosen for classes of diagrams. In addition, Monte Carlo integration of the loops and over the phase-space of real emissions can be merged, potentially stabilising and speeding up the calculation.

Appendix B: The one-point FDR integral
We compute the one-loop integral By taking the μ → 0 limit at the integrand level and replacing μ with μ R in the logarithmic divergent vacuum one obtains (B.14) with R 0 and R ν defined in (123). One computes .
Note that there is no pole in this case. In Table 13 we report a comparison between a numerical implementation of (B.16) and the analytic result A FDR /m 2 = iπ 2 (1 + ln ν). (B.17)