Monte Carlo evaluation of divergent one-loop integrals without contour deformation

Reference [1] introduces a method for computing numerically four-dimensional multi-loop integrals without performing an explicit analytic contour deformation around threshold singularities. In this paper, we extend such a technique to massless scalar one-loop integrals regularized in the framework of dimensional regularization. A two-loop example is also discussed.


Introduction
In recent years a huge effort has been devoted to the problem of the computation of loop integrals.The reason is the need of accurate theoretical predictions able to cope with the ever-increasing precision of the data collected in particle physics experiments.
Regardless of the approach employed for regularizing infinities, the finite part of the calculation is plagued by the presence of the so-called threshold singularities.These are integrable singularities avoided by the +iϵ a e-mail: pittau@ugr.espropagator prescription.Such structures are not a problem for analytic calculations, but must be properly addressed when using numerical techniques.
In a previous paper [1], a method has been introduced, which permits an accurate fully numerical treatment of threshold singularities that can be easily implemented in Monte Carlo (MC) codes.The advantage of this technique is that, even if a non-zero ϵ must be kept, its influence on the result can be lowered close to the machine precision level, e.g. between 10 −12 and 10 −9 times the largest physical scale appearing in the problem.
The performance of this procedure has been studied in [1] in the case of finite multi-loop Feynman integrals, or divergent ones regularized via four-dimensional methods.In this paper we extend for the first time this approach to integrals regularized within dimensional regularization, focusing our attention on scalar integrals that provide a complete basis for any one-loop calculation in massless theories [36,37].Higher-loop integrals will be studied in detail elsewhere, although we envisage that the experience gathered at the one-loop level can be comfortably adapted to multi-loop environments.
The structure of the paper is as follows.In Sect. 2 we briefly review the method.Section 3 fixes our kinematics and conventions.Sections 4, 5 and 6 are devoted to the study of the massless 2-, 3-and 4-point scalar oneloop integrals, respectively, and in Sect.7 we present a simple two-loop example.
Finally, it is important to mention that throughout the paper we distinguish between the ϵ and ε symbols.The latter parameterizes the n-dimensional loop integration, n = 4 − 2ε, while the former denotes the contour deformation around single-pole singularities.All results presented in this paper are produced with ϵ = 10 −9 .

The method
In this section we briefly recall the method of [1].Our aim is to flatten the singular behavior of a threshold singularity parameterized as where the numerator function F (x) is regular in x = 0. To achieve this we introduce a complex integration variable z = α + iβ related to x by x + iϵ = e iπ (1−z) .
The requirement that x remains real fixes the path in the complex z plane to be so that Using now gives Two comments are in order.Firstly, the integrand of ( 5) is now regular in x = 0 for arbitrarily small values of ϵ.In fact, the ϵ dependence is moved to the boundaries of the integration region, x = ±1, which are reached exactly only when ϵ → 0. However, x = ±1 are far away from the threshold singularity of (1).That explains why the algorithm survives tiny numerical values of ϵ.Secondly, if F (x) contains branch cuts in the x complex plane, the fact that x always lies on the real axis ensures that the right Riemann sheet is automatically taken when −1 ≤ x ≤ 1.Thus, compared to methods based on contour deformation, one does not have to worry about choosing a path that avoids the pole at x = −iϵ without crossing any cut of F (x).
Equation ( 5) is optimal for integrating over α.To flatten the integral over β, the parameterization complementary to (4) is needed, namely However, (2) implies that α is a two-valued function of β.Therefore, it is necessary to divide (5) into two parts where y α := ϵ/tan(απ).Inserting ( 6) into (7) gives with that is optimized for the integration over β.In practice, equations ( 7) and ( 8) can be merged together by means of a multichannel MC approach [38],1 so that the complete 1/(x + iϵ) behavior of ( 1) is flattened.The described algorithm is implemented in the code GLoop [39].The present version is able to deal with integrals of the type with m up to 4. The numerical results presented in this paper require m = 3 at most.

Kinematics and loop integration
For the purposes of this work it is sufficient to consider a p 1 + p 2 → p 3 + p 4 massless kinematics given by In (10) s = (p 1 + p 2 ) 2 and θ 13 is the scattering angle defined by the relation The n-dimensional loop momentum is where c θ = cos θ, s θ = sin θ, c ϕ = cos ϕ, with 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ 2π.
Rescaling p 1,2,3 and q by √ s produces the following dimensionless vectors span a 3-dimensional space, so that in the one-loop integrands one can trade ω α n for its 4-dimensional projection defined as In fact The scalar 2-point one-loop function of (20).
Finally, the integration over the rescaled loop momentum can be parameterized as where ´n is the (n − 1)-dimensional integration volume.
In terms of ρ, θ and ϕ it reads ˆn = 2π If the integrand is independent of ϕ, the integration over c ϕ can be carried out analytically by using that gives ˆn = 2π Likewise, integrating over c θ produces ˆn = 2π 4 The UV divergent one-loop 2-point integral As a first illustration of our procedure we compute by MC the dimensionally regularized one-loop scalar integral of Fig. 1, where D 0 = q 2 + iϵ and D 1 = (q + p 1 + p 2 ) 2 + iϵ.Given that B(s) diverges in four dimensions, our strategy is to split it into a finite part and a UV divergent piece, in a way that the former can be computed numerically using the 4-dimensional algorithm of Sect.2, while the latter is evaluated analytically.
Rescaling loop and external momenta by √ s gives with d 0 = τ 2 − ρ 2 + iϵ and d 1 = (τ + 1) 2 − ρ 2 + iϵ.The integrand of (22) does not depend on θ and ϕ.Hence, ´n can be taken as in (19).We now integrate over τ by using the residue theorem.The result is To achieve the splitting of ( 21), we subtract and add back an integrand with the same ρ → ∞ behavior of (23).Among the various possibilities we choose 1 ρ(ρ 2 +1/4) .This allows us to recast B = B F + B UV , in which the finite piece reads while B UV is easily computed analytically, Introducing the variable σ = ρ 2 − 1/4 allows one to rewrite (24) in a form suitable to be integrated with GLoop, where Θ is the Heaviside step function.
Our MC estimate with 10 7 MC shots gives B F /(iπ 2 ) = 3(9) × 10 −4 + i 3.1414 (7), (27) to be compared to the analytic value The time to produce the result of ( 27) on a single 2.2 GHz processor is of about 1.6 s.
Finally, the analytic continuation to s < 0 is achieved by replacing s/µ 2 → s/µ 2 + iϵ in (22) [40].This produces with Lastly, it is noteworthy to mention that from a numerical standpoint, the treatment of UV divergences is remarkably similar in terms of dimensional regularization and FDR [27][28][29][30][31]33].The sole distinction lies in their action towards the subtracted UV divergent part, which is computed analytically and added back in the former approach, whereas it is judiciously discarded in FDR [1].
5 The IR divergent one-loop 3-point integral Here we study the dimensionally regularized triangle function of Fig. 2 with two massless momenta p 2 1,2 = 0 and massless denominators As in the previous section, we rescale all dimensionful quantities by √ s.This produces The rescaled integral C reads with ´n given in (18) and C is divergent in the soft and collinear limits, namely it develops 1/ε and 1/ε 2 poles under the n-dimensional integration.To be able to perform the loop integration numerically, we first construct an approximation C IR whose integrand subtracts the infrared behavior in a local fashion.Then we reinsert the result of an analytic computation of C IR .Schematically, C = C F + C IR , where is computed by MC.
Using the residue theorem to integrate over τ gives Only the first two terms of ( 35) are divergent when ρ < 1/2.An approximation sharing their IR behavior is constructed by expanding R around ρ(1 + c θ ) = 0, that produces which can be easily evaluated analytically The integrand of ( 36) can now be subtracted to (35) and the resulting integral produces the finite contribution of (34).In four dimensions (18) gives with R ± = 1/2 ± ρ, thus where the −iϵ is kept only in denominators with threshold singularities.Now we put (38) in a form suitable to be integrated with GLoop.The terms between curly brackets have denominators of the form R + r i , where the r 1÷7 are independent of R. We then introduce two integration variables defined as follows and rewrite The numerator F C (σ 1 , σ 2 ) of ( 40) is fully expressible in terms of Heaviside functions and is given in Appendix A. With 4 × 10 8 MC points our estimate is to be compared to the analytic result The time to produce 10 6 MC shots on a single 2.2 GHz processor is of about 0.32 s.In our computation we have assumed s > 0. The analytic continuation to s < 0 is again obtained by replacing s → s + iϵ in (31).Expanding in ε gives where L = ln − s/µ 2 − iϵ .
6 The IR divergent one-loop 4-point integral In this section we consider the massless box diagram of Fig. 3, where p 2 1,2,3,4 = 0.By rescaling all dimensionful quantities by √ s one arrives at where we have defined x = −t/s, so that in the physical region one has 0 ≤ x ≤ 1.The rescaled 4-point integral reads with ´n given in (16).The denominators d 0,1,2 are as in (33).Furthermore, As before, to compute (46) by MC, we first have to subtract a simpler function D IR (x) with the same local IR behavior of D(x).We choose that can be integrated analytically by means of ( 42) and ( 43), The rationale behind (48) is as follows.On the one hand, it is well known that the four 3-point integrals produce the same 1/ε and 1/ε 2 poles of D(x) [41].On the other hand, the last term is chosen in such a way that it compensates their finite contribution.In this way gives the finite part of D(x) directly.
Integrating over τ the IR finite combination of integrals appearing in (49) gives where P denotes the Cauchy principal value, To derive (50) we have systematically identified terms related by the interchange ρ ↔ S.This is possible because ´4 is invariant under shifts and rotations of the spatial components of the vectors π α 1,2,3 and ω α .The ϕ dependence is entirely contained in S 2 − ρ 2 and can be integrated out by using Finally, inserting (52) and ( 50) into (49) gives with Our MC estimates are presented in Table 1 and compared to the analytic result [40] The time to generate 10 6 MC shots on a single 2.2 GHz processor is of about 0.4 s.

A two-loop example
In this section we compute the two-loop bubble-box diagram of Fig. 4, which has collinear and soft IR divergences in addition to a UV-divergent sub-diagram, . More precisely, we evaluate its IR and UV subtracted counterpart obtained as described and reported in [42], The first term of (56) is the original integral, while are the fractions of the momenta p 1 and p 2 carried by the internal lines with momenta q and p 1 + p 2 − q, respectively.Finally, m is an arbitrary mass used to subtract the UV behavior.
The best way to apply the algorithm of Sect. 2 to the case at hand is to use the gluing procedure described in [1], which allows one to express (56) in terms of a treelevel part convoluted with the one-loop sub-diagram.To achieve this, we perform analytically the k integration and rescale all dimensionful quantities by √ s, which produces where µ 0 = m 2 /s and In (58) d 0,1,2 = D 0,1,2 /s are rescaled denominators, Next, we trade the integrations over τ , ρ and c θ for integrations over σ 0,1,2 defined in (58).This gives with where λ(x, y, z) is the Källén function.The integration over the azimuthal angle ϕ can be performed by using ˆ2π After that, T F (x) can be written in terms of a triple integral suitable to be evaluated numerically with GLoop, with F T (σ 0 , σ 1 , σ 2 , x) given in Appendix C. Note that the numerator F T contains branch cuts controlled by the iϵ prescription (see (62)), but because our algorithm maintains σ 0,1,2 in the real axis, the correct Riemann sheet is automatically taken.The analytic result in the physical region 0 ≤ x ≤ 1 reads [42] T where S 12 is the Nielsen polylogarithm.Fig. 4 The two-loop bubble-box of (55).

Conclusion and outlook
Any attempts towards a numerical loop integration requires controlling threshold singularities.A possible approach is contour deformation [43,44], that calls for analytic knowledge of the cut structure of the integrand [45] or numerical checks establishing whether the deformation stays on the correct side of the singularity.In [1] an alternative has been proposed and shown to be effective in the MC estimate of four-dimensional multi-loop integrals directly in Minkowski space.
In this paper we have extended this technique to massless scalar one-loop integrals with no more that four external legs, regularized within dimensional regularization.Our strategy is based on a separation of the 1/ε and 1/ε 2 poles before integration.The finite part, containing the threshold singularities, can then be integrated numerically in four dimensions.A fully numerical evaluation of one-loop amplitudes with more than four legs is in principle possible if the coefficients of the contributing one-, two-, three-and four-point integrals are also determined in a numerical fashion by using, for instance, the method of [37].
We have presented numerical results obtained with the help of the code GLoop.A MC error of the order of a few per mil can usually be obtained for a modest CPU cost.As with any other numerical method, this level of precision is expected to be sufficient for phenomenological purposes when the gauge cancellations are moderate.If this is not the case, cancellations among diagrams must be enforced to occur before integration.This should be feasible because they are usually controlled by Ward identities operating at the integrand level.
Enlarging the range of applicability of our method beyond one loop requires removing UV and IR divergences by adding appropriate counterterms at the level of the integrand.Ultraviolet counterterms can be constructed by using, for instance, the same procedure that defines FDR integrals [27][28][29][30][31]33].As for the infrared behavior, a systematic approach has been developed for two-loop integrals by Anastasiou and Sterman in [42].In Sect.7 we have presented a simple two-loop example showing how the subtraction method of [42] can be combined with our algorithm.For all these reasons, we believe that the strategy described in this paper can be extended up to two loops.A deeper exploration of this issue is planned for a future publication.
Table 2 shows a comparison between our MC estimate based on (63) and (64).The time to produce 10 6 MC shots on a single 2.2 GHz processor is of about 0.42 s.

Table 2
Numerical estimates of T F (x) in (63) for several values of x = −t/s and µ 0 = 1.The analytic result is reported in (64).Numbers obtained with 8×10 10 MC shots.MC errors between parentheses.