Direct numerical evaluation of multi-loop integrals without contour deformation
Abstract
The work of RP is supported by the SRA grant PID2019-106087GB-C21 (10.13039/501100011033), by the Junta de Andalucía grants A-FQM-467-UGR18 and P18-FR-4314 (FEDER), and by the COST Action CA16201 PARTICLEFACE. The work of BW was partially supported by STFC HEP consolidated grants ST/P000681/1 and ST/T000694/1.
Full text
Eur. Phys. J. C (2022) 82:55 https://doi.org/10.1140/epjc/s10052-022-10008-6 Regular Article - Theoretical Physics Direct numerical evaluation of multi-loop integrals without contour deformation Roberto Pittau1,2,a, Bryan Webber3,b 1Departamento de Física Teórica y del Cosmos and CAFPE, Universidad de Granada, 18071 Granada, Spain 2Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland 3University of Cambridge, Cavendish Laboratory, J.J. Thomson Avenue, Cambridge, UK Received: 1 November 2021 / Accepted: 6 January 2022 © The Author(s) 2022 Abstract We propose a method for computing numerically integrals defined via ideformations 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 fourdimensional regularization techniques, our treatment can be extended to ultraviolet and infrared divergent integrals. 1 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, ae-mail: [email protected] (corresponding author) be-mail: [email protected].cam.ac.uk where single poles are moved away from the real integration domain by the iprescription. These singularities require special treatment, such as a contour deformation into the complex plane [2–5], or vanishing-width extrapolations methods [6–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 iprescription 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–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 2details our approach. In Sect. 3we use it to integrate numerically threshold singularities after an analytic integration over the energycomponents of the loop momenta. Section 4explains how to glue together lower-loop structures to compute numerically certain classes of higher-loop MIs. Finally, in Sects. 5and 6 we extend our treatment to ultraviolet and infrared divergent configurations regularized via the four-dimensional method of [13]. 0123456789().: V,-vol 123
55 Page 2 of 22 Eur. Phys. J. C (2022) 82:55 2 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. 2.1 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 xwhen −1≤x≤ 1, yif 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 I:= lim →0ˆ1 −1 dx 1 x+i(1) is flattened with x∈R. This is obtained by imposing x+i=eiπ(1−z),(2) where zis a new complex integration variable. In fact, inserting (2)in(1) gives the desired result, I=−iπ$1 0 dz.(3) Equation (3) evaluates to −iπalong any curve in the zcomplex plane connecting z=0toz=1 when →0. We use this freedom to impose x∈Rby parametrizing z=α+iβwith α, β ∈Rand π≤α≤1− π.(4) Inserting (4)in(2)gives x=eπβ cos[π(1−α)]+i{eπβ sin[π(1−α)]−},(5) which is real when πβ =ln sin[π(1−α)], namely x=xα:= tan[π(1−α)].(6) Therefore dz =dα1+idβ dα,(7) which gives $1 0 dz =lim →0 1 gˆ1−/π /π dα1+ixα ,g:= 1−2 π, (8) where ghas 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)toanintegration over a function f(x)=φ(x)/(x+i), (9) with φ(x)sufficiently smooth at x=0,1 If:= ˆ1 −1 dx f(x)=ˆ1 0 dyf(−y)+f(y).(10) Splitting the integration region of (8) into the two sectors with xα<0orxα>0gives If=−iπ gˆ1/2 /π dα ×1−iyα φ(−yα)+1+iyα φ(yα),(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), If=ˆ1 0 dρf(−y)+f(y) g(y),(12) with ´1 0dy g(y)=1. By comparing (12)to(11) one determines g(y)=2 π(y2+2),y= tan(απ),α= π+ρg 2.(13) 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, dz =dβdα dβ+i,(14) which gives If=−iπ gˆβ+ β− dβ × −yβ+iφ(−yβ)− yβ+iφ(yβ),(15) where yβ:= eπβ1− eπβ 2 ,β −=1 πln sin ,β +=ln π. 1From now on, we omit lim→0and consider as an infinitesimal parameter. 123
Eur. Phys. J. C (2022) 82:55 Page 3 of 22 55 Again, (15) is correctly normalized also for a small but not vanishing . Comparing (12)to(15) gives now g(y)=− g ln(sin ) y (y2+2),y=eπβ1− eπβ 2 , β=ln() −ρln(sin ) π.(16) Multichanneling Flattening the whole 1/(x+i) behaviour of (9) requires a merging of (13) and (16), whose densities we dub g1(y) and g2(y). This can be achieved via a multichannel approach with combined density gc(y):= α1g1(y)+α2g2(y)and α1+α2=1, If=ˆ1 0 dρf(−y)+f(y) gc(y).(17) In (17), ρis generated according to the distribution g1,2(y) with probability α1,2, and the a-priori weights α1,2can be optimized as described in [14]. 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 gi(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 ˆ1 0 dygi(−y)+gi(y)=ˆ1 −1 dx gi(x), (18) which means that if xis randomly chosen in −1≤x≤1, the density is gi(−|x|)+gi(|x|). Hence, the MC weight is f(−y)+f(y)/gi(−y)+gi(y), with ginormalized such that ˆ1 −1 dx gi(x)=1.(19) In summary, with Nch channels (including g1and g2), the more general multichannel MC mapping reads If=ˆ1 0 dρf(−y)+f(y) gtot(y),dρ=gtot(y)dy,(20) where gtot(y)=gc(y)+ Nch i=3 αi(gi(−y)+gi(y)), (21) with arbitrary (but self-adjustable) weights fulfilling Nch i=1αi=1. In the actual MC used to produce the results presented in this paper we superimpose on gc(y)a flat distribution g3(x)=1/2 and a channel g4(x)=1 2ln 1+δ δ 1 1−|x|+δ,δ=10−4,(22) 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 − ˆdσ:= lim →∞ˆ − dσ, (23) which can be mapped onto the interval [−1,1]by changing variable, σ=x 1−x2.(24) Thus, for instance, − ˆdσφ(σ) σ+i=ˆ1 −1 dx 1+x2 1−x2φx 1−x21 x+i,(25) 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 If,n:= ˆ1 −1 n j=1dxjf({x}), (26) with f({x})=φ({x})/n j=1(xj+i). Our notation is such that {x}=x1,x2,...,xnand φ({x})is a smooth function at {x}={0}. The result is If,n=ˆ1 0 n j=1dρjf(−{y})+f({y}) n j=1gtot(yj),(27) where dρj=gtot(yj)dyjand the numerator stands for a sum over the 2nterms with positive or negative arguments. For instance, when {y}=y1,y2, f(−{y})+f({y}) =f(−y1,−y2)+f(−y1,y2) +f(y1,−y2)+f(y1,y2). (28) 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. 123
55 Page 4 of 22 Eur. Phys. J. C (2022) 82:55 2.2 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 I1[φ]:=ˆa −a dx φ(x) x+i=ˆa 0 dx xφ−(x)−iφ+(x) x2+2, (29) where φ±(x)=φ(x)±φ(−x). We can write this as I1[φ]=ˆa/ 0 dy yφ−(y)−iφ+(y) 1+y2 =ˆrm 0 dr φ−(e2r−1)−iˆθm 0 dθφ +( tan θ), (30) where rm:= ln(1+a2/2)/2,θ m:= arctan(a/). (31) 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 /(x2+2)and x/(x2+2)behaviours of (29) correspond to the local densities in (13) and (16), respectively. The method is easily generalised to two variables. Consider I2[φ]:=ˆa −a dx1dx2 φ(x1,x2) (x1+i)(x2+i) =ˆa 0 dx1dx2 (x1,x2;) (x2 1+2)(x2 2+2),(32) where (x1,x2;) =x1x2φ11 −i(x1φ10 +x2φ01)−2φ00 (33) with φ00 =φ(x1,x2)+φ(−x1,x2)+φ(x1,−x2)+φ(−x1,−x2), φ10 =φ(x1,x2)−φ(−x1,x2)+φ(x1,−x2)−φ(−x1,−x2), φ01 =φ(x1,x2)+φ(−x1,x2)−φ(x1,−x2)−φ(−x1,−x2), φ11 =φ(x1,x2)−φ(−x1,x2)−φ(x1,−x2)+φ(−x1,−x2). (34) For MC evaluation, we proceed as follows: for each shot, generate x1r,x1t,x2r,x2twhere x1r=e2r1−1,x1t=tan θ1, x2r=e2r2−1,x2t=tan θ2(35) where 0<r1,2<rm,0<θ 1,2<θ m(36) uniformly, with rmand θmas in (31). In (34), set x1=x1t when the first superscript is 0 and x1=x1rwhen it is 1, and similarly for x2according to the second superscript. The weights for the real and imaginary parts are then wr=φ11r2 m−φ00θ2 m,w i=−(φ10 +φ01)rmθm.(37) The generalisation of (32)tonvariables is clear: φ{kj}has superscript kj=1inthe jth location when there is an xjin the integrand, otherwise kj=0. The symmetrized function becomes ({xj};) =(−i)n {kj=0,1} n j=1 (ixj/)kjφ{kj}({xj}), (38) where φ{kj}({xj})= {lj=0,1} n j=1 (−1)kjljφ{(−1)ljxj}.(39) For MC evaluation, for each shot, generate two points in the n-dimensional hypercube xjr =e2rj−1,xjt =tan θj,(40) where again 0 <rj<rmand 0 <θ j<θ muniformly. In (39), set xj=xjt when kj=0 and xj=xjr when kj=1. The weight is then wr+iwi=(−iθm)n {kj=0,1} n j=1 (irm/θm)kjφ{kj}({xj}). (41) Note that each shot involves 2nrandom numbers for {xjt} and {xjr}, and then 4nfunction evaluations at xj=±xjt and ±xjr, so the computation time increases rapidly with the number of variables. 2.3 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 T() =ˆ1 −1 3 j=1dxj xj+i1 j1,j2,j3=0 xj1 1xj2 2xj3 3(42) =(8−6π2)+iπ(π2−12), (43) whose behaviour at xj∼0 mimics a typical multidimensional environment. The result of this comparison is given in Fig. 1, where the solid (dashed) line represents the real (imaginary) part of (43). Bullets and squares with errors are the MC predictions for e[T()]and m[T()], respectively. To quantify the effect of a nonzero on the MC esti123
Eur. Phys. J. C (2022) 82:55 Page 5 of 22 55 Fig. 1 MC results for the real (red bullets) and imaginary (blue squares) part of (42). They have been obtained with 1010 MC shots per point, corresponding to 8×1010 calls to the integrand. To minimize the statistical fluctuations, the same sequence of random numbers is used for all values of mate QMC() ±ΔQ() of a known quantity Q,itisconvenient to introduce the estimators Δ1() =|Q−QMC()| |Q|,Δ 2() =max , ΔQ() |Q|.(44) Requiring the = 0biasonQMC() to be of the order of the maximum between and the relative MC error gives the condition R(, Q):= Δ1() Δ2() ∼O(1). (45) Table 1reports the Rvalue of the entries of Fig. 1and their MC accuracy defined as δ() := max Δe[T()],Δm[T()] |TMC()|.(46) From Fig. 1and Table 1we infer that a range 10−8≤≤ 10−6is adequate to achieve MC estimates accurate at the level of three parts in 105. 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 1shows 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. 3 A semi-numerical integration algorithm for MIs In this section we illustrate how the approach of Sect. 2can be successfully applied to produce stable and precise seminumerical MC estimates of loop MIs. This is achieved in Table 1 The ratio in (45) for the real and imaginary parts of T() in (42) computed with 1010 (3×1011) MC shots when ≥10−8(=10−12). The last column reports the MC accuracy defined in (46) R(, e[T])R(, m[T]) δ() 10−41.1 7.0 3×10−5 10−51.1 0.5 3×10−5 10−60.6 0.3 3×10−5 10−70.5 0.2 3×10−5 10−80.6 0.1 3×10−5 10−12 0.02 1.5 7×10−6 Fig. 2 The scalar three-point one-loop function with arbitrary kinematics and masses C(P2,p2 1,p2 2,m0,m1,m2) 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 iprescription 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.2In what follows we illustrate the performance of this strategy by means of two examples. 2The use of threshold counterterms [15] could further improve the precision of our approach. 123
55 Page 6 of 22 Eur. Phys. J. C (2022) 82:55 Table 2 Numerical estimates of the one-loop integral (52) multiplied by m2, compared to the analytic result of [16]. Numbers obtained with 2 ×107 MC points. The MC errors are indicated between parentheses √τMC result Analytic result 0.01 9.85(2)×10−7−i4.9408(92) 0 −i4.9348 0.2 9.90(1)×10−7−i4.9482(55) 0 −i4.9513 0.5 1.0233(8)×10−6−i5.0361(41) 0 −i5.0412 1.99 1.4341(9)×10−5−i1.0783(3)×1010−i1.0782×101 2.01 1.5350(6) −i1.2006(3)×1011.5343 −i1.2006×101 10 1.4216(3) +i5.5007(30)×10−11.4216 +i5.5030×10−1 1022.8562(8)×10−2+i3.6999(15)×10−22.8557×10−2+i3.6990×10−2 1045.7141(20)×10−6+i1.6258(5)×10−55.7116×10−6+i1.6258×10−5 3.1 A one-loop example Consider the three-point function of Fig. 2in the case m0= m1=m2=m,p2 1=p2 2=0 and timelike P. Rescaling all momenta by m q m=(t,ρcθ,ρsθsφ,ρsθcφ), P m=(√τ,0,0,0), (47) p1 m=√τ 2(1,1,0,0), p2 m=√τ 2(1,−1,0,0), (48) gives C0:= C(P2,0,0,m,m,m)=2π m2ˆ1 −1 dcθˆ∞ 0 dρρ2 ׈∞ −∞ dt 1 (σ0+i)(σ1+i)(σ2+i),(49) with σ0=q2 m2−1,σ 1=(q−P)2 m2−1,σ 2=(q−p2)2 m2−1. One splits 1 σ2+i=1 2R21 t−√τ/2−R2+i −1 t−√τ/2+R2−i(50) with R2 2:= ρ2+τ 4+√τρcθ+1. Thus ˆ1 −1 dcθ 1 σ2+i=1 √τρln t−√τ/2−R− 2+i t−√τ/2−R+ 2+i +ln t−√τ/2+R− 2−i t−√τ/2+R+ 2−i,(51) where R± 2:= (√τ/2±ρ)2+1. The cut of the logarithms with +i(−i) is in the lower (upper) tcomplex half-plane, so that the integration over tin (49) is trivial once one rewrites 1 (σ0+i)(σ1+i) =1 4R2 01 t−R0+i−1 t+R0−i Fig. 3 The two-loop self-energy diagram S2(m,m0,m1) ×1 t−√τ−R0+i−1 t−√τ+R0−i, with R2 0:= ρ2+1. The results is C0=2iπ2 m2τˆ∞ −∞ dr r−iθr+√τ 2−1Lr+√τ 2,√τ +θ(r−√τ 2−1)L(r−√τ 2,−√τ),(52) where L(r,√τ) := ln r−√τ 2+R(r,−√τ)−i r−√τ 2+R(r,√τ)−i and R(r,√τ) := τ/4+r2+√τr2−1. When √τ>2, the first integrand of (52) develops a pole at r=ithat migrates towards the integration region in the limit →0. Treating this with the strategy of Sect. 2gives the results presented in Table 2. 3.2 A two-loop example We study the two-loop self-energy scalar diagram of Fig. 3. For a timelike P/m=(√τ,0)and m1=mit reads S2:= S2(m,m0,m)=1 m2ˆd4ω1d4ω2 5 j=1 1 σj+i,(53) 123
Eur. Phys. J. C (2022) 82:55 Page 7 of 22 55 where ωi:= qi/m=(ti,ρi)(54) and σ1=ω2 1−1,σ 2=ω2 2−1, σ3=σ1+τ−2√τt1,σ 4=σ2+τ+2√τt2, σ5=(t1+t2)2−ρ2 1−ρ2 2−2ρ1ρ2cθ−μ0,(55) with μ0:= m2 0/m2. Integrating over the angular variables gives S2=4π2 m2ˆ∞ 0 dρ1ρ1ˆ∞ 0 dρ2ρ2ˆdt1dt2 4 j=1 1 σj+i ×ln (t1+t2)2−(ρ1−ρ2)2−μ0+i (t1+t2)2−(ρ1+ρ2)2−μ0+i.(56) The integration over t1and t2is trivial and produces S2=2π4 m2τ λ1,2=±ˆ∞ −∞ dr1ˆ∞ −∞ dr2 F(r1,r2,λ 1,λ 2) (r1−i)(r2−i),(57) where F(r1,r2,λ 1,λ 2)=λ1λ2θ(A1−1)θ(A2−1) ×ln r1+r2+(A2 1−1+A2 2−1)2+μ0−i r1+r2+(A2 1−1−A2 2−1)2+μ0−i and Ai:=ri+λi√τ/2. Threshold singularities at r1,2=i are present when λ1,2=+1if√τ>2. When m0=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 109MC-point estimate with ρ=.1gives −m2τ π4S2(m,m,m)=8.582(6)−i2.706(4). (58) 4 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, ω:= q/m=(t,ρcθ,ρsθsφ,ρsθcφ). (59) Furthermore, we define μi:= m2 i/m2,τ:= P2/m2,χ := (p2−p3)2/m2, τi:= p2 i/m2,τ ij := τi−τj, λij := λ(τ, τi,τj), Table 3 The two-loop integral (57) with m0=0 multiplied by −m2τ/π4for several values of ρ:= 4/τ. Numbers obtained with 109(1010) MC points when ρ>1(ρ<1). The analytic result is taken from [17]. MC errors between parentheses ρMC result Analytic result 0.1 8.49(1) −i1.94(2) 8.495 −i1.927 0.3 9.34(1) −i5.47(2) 9.340 −i5.460 0.5 9.19(1) −i9.71(1) 9.195 −i9.716 0.7 7.39(1) −i15.79(1) 7.396 −i15.783 0.9 −1.03(2) −i27.591(8) −1.061 −i27.581 1.1 −15.538(2) −i1.8314(4)×10−5−15.540 +i0 1.3 −7.9915(8) −i5.1218(7)×10−6−7.9921 +i0 1.5 −5.5608(6) −i2.9000(4)×10−6−5.5614 +i0 1.7 −4.2990(5) −i2.0139(3)×10−6−4.2996 +i0 1.9 −3.5153(5) −i1.5412(2)×10−6−3.5157 +i0 k2:= λ12/τ2,k±:= k2±i, (k)2:= λ34/τ2, and study cases up to a P:= p1+p2→p3+p4kinematics of the form p1=m 2√τ(τ +τ12,λ 1 2 12,0,0), p2=m 2√τ(τ −τ12,−λ 1 2 12,0,0), p3=m 2√τ(τ +τ34,λ 1 2 34 cos θ13,λ 1 2 34 sin θ13,0), p4=m 2√τ(τ −τ34,−λ 1 2 34 cos θ13,−λ 1 2 34 sin θ13,0). (60) Rescaled propagators belonging to the loop momentum qare denoted by σ0:= q2/m2−μ0,σ 1:= (q−P)2/m2−μ1, σ2:= (q−p2)2/m2−μ2,σ 3:= (q−p3)2/m2−μ3. In addition, we define λ:= λ(τ, σ0+μ0,σ 1+μ1), (61) and3 σij a:= τ1−τ+σj+μj+kcθ 2λ1 2(τ, σi+μi,σj+μj) +1 21−τ12/τ(τ +σi+μi−σj−μj), σij b:= σi+μi+σj+μj+τ1+τ2−τ−σij a, σij c:= τ3+σi+μi+k 2λ1 2(τ, σi+μi,σj+μj) 3σij a,b,c,dare the invariants (q−p1,2,3,4)2/m2computed at values of t and ρsatisfying the conditions t2−ρ2=σi+μiand (t−√τ)2−ρ2= σj+μj. 123
55 Page 8 of 22 Eur. Phys. J. C (2022) 82:55 ×cθcos θ13 +sθsφsin θ13 −1 21+τ34/τ(τ +σi+μi−σj−μj), σij d:= σi+μi+σj+μj+τ3+τ4−τ−σij c.(62) The essence of the procedure is to use σ0and σ1as integration variables of the method of Sect. 2. This is achieved by multiplying the integrand by 1=− ˆdσ0− ˆdσ1Δ(σ0,σ 1,ρ,t), (63) where Δ(σ0,σ 1,ρ,t):= δ(σ0+μ0+ρ2−t2) ×δ(σ1−σ0+μ1−μ0−τ+2√τt). (64) This gives rise to the appearance of the following three functionals, (μ0,μ1) 1[J1]:=ˆd4ωJ1Δ(σ0,σ 1,ρ,t), (μ0,μ1,μj) j[Jj]:=ˆd4ωJj σj+iΔ(σ0,σ 1,ρ,t), (65) where j=2,3. Assuming J1independent of any angular variable and J2(J3) independent of θ(φ) allows one to compute the functionals once for all. 1reads (μ0,μ1) 1[J1]= π 2τλ1 2θ(λ)J1.(66) As for 2, one has (μ0,μ1μ2) 2[J2]= 1 4τθ(λ)1 k+ ln 1+k+λ1 2 A2+i −1 k− ln 1−k−λ1 2 A2+iˆ2π 0 dφJ2. (67) with A2:= (σ0+μ0)1+τ12 τ+(σ1+μ1)1−τ12 τ +τ1+τ2−τ−2μ2.(68) Note that (66) and (67) have been analytically continued to configurations with any sign of τ,τ1,τ2,λ12. Finally (μ0,μ1,μ3) 3[J3]= π 2τλ1 2θ(λ)ˆ1 −1 dcθ J3 A3 1 1−B3 A2 3 ,(69) in which A3:= (σ0+μ0)1−τ34 τ+(σ1+μ1)1+τ34 τ +τ3+τ4−τ−2μ3+λ1 2kcθcos θ13 +i, B3:= λλ34 τ2s2 θsin2θ13.(70) Fig. 4 The scalar four-point one-loop function D(p2 1,p2 2,p2 3,p2 4,(p1 +p2)2,(p2 −p3)2,m2 1,m2 2,m2 0,m2 3)with arbitrary kinematics and masses 4.1 One-loop examples To elucidate the procedure, we first consider gluing tree-level structures to compute the three-point function of Fig. 2–with arbitrary kinematics and masses–and the box diagram of Fig. 4with p2 i=0 and m0=m1=m2=m3=m, which we dub Cand D0, respectively. Equations (63), (67) and (69) produce m2C=− ˆ1 j=0dσj σj+i(μ0,μ1) 2[1],(71) m4D0=− ˆ1 j=0dσj σj+i(μ0,μ1,μ3) 31 σ2+i.(72) In (71) the integration over the azimuth angle φis trivial, (μ0,μ1) 2[1]= π 2τθ(λ) j=±1 jkj ln 1+jkjλ1 2 A2+i,(73) while the integral over cθin (72) has to be dealt with numerically using the method of Sect. 2, which gives (μ0,μ1,μ3) 31 σ2+i=π τθ(λ0)− ˆdσ2 σ2+i ×θ(σ2−σ− 2)θ(σ + 2−σ2) A01−B0 A2 0 ,(74) where σ± 2=στ±λ 1 2 0 2,σ τ:= σ0+σ1−τ, A0=2σ2−χ τ(στ−2σ2)+i, B0=−16χ τ1+χ τ[σ2(στ−σ2)−τ−σ0σ1], λ0=λ(τ, σ0+1,σ 1+1). (75) 123
Eur. Phys. J. C (2022) 82:55 Page 9 of 22 55 Improving the numerical accuracy When inserted in (71) and (72), Eqs. (73) and (74) could potentially produce inaccurate results when strong cancellations are expected among different integration regions. This happens if (a) the integrands do not vanish fast enough at large values of σ0,1; (b) τis small. (76) Note that case (a) is relevant to Cbut not to D0, since lim σ0,1→∞(μ0,μ1) 2[1]∼constant ,(77) lim σ0,1→∞(μ0,μ1,μ3) 31 σ2+i∼1 σ0,1 ,(78) while (b) applies to both Cand D0due 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 Cfinite despite (77).4We define σ10 := σ1−σ0,(79) in terms of which m2C=lim →∞ˆ − dσ10 ˆdσ0 (μ0,μ1) 2[1] (σ0+i)(σ0+σ10 +i). Now the σ0integral is convergent by power counting and λ in (61) behaves as lim σ10→∞λ∼˜ λ:= (σ10 +β)2−α, (80) where αand βare constants. Replacing λwith ˜ λin (73) produces an integrand in which all branch points and poles are located in the lower σ0complex half-plane. As a result, the integral over σ0approaches zero when σ10 →±∞,so that the →∞limit exists. This same reasoning allows one to construct a class of vanishing integrals defined as m2˜ C(α,β, 0):= ˆd4ωθ(|σ10|−0) (σ0+i)(˜σ1+i)(˜σ2+i),(81) with m2˜σ1:= q2+P2−m2 1−2(˜q·P), m2˜σ2:= q2+p2 2−m2 2−2(˜q·p2), (82) 4The principal value integrations in (71) are not sufficient to regularize the large σ0,1behaviour. In fact, (μ0,μ1) 2[1]approaches two different constants when σ0,1→∞or σ0,1→−∞, so that no cancellation is possible. where ˜qis the asymptotic |σ10|→∞limit of q, ˜q m:= (˜ t, ˜ρ), ˜ t:= −σ10 +β 2√τ,˜ρ:= ˜ λ1 2 2√τ.(83) Again, all cuts and poles lie in the lower σ0complex halfplane, so that ˜ C(α,β, 0)=0.(84) Now αand βcan be set to obtain a local cancellation of the problematic large σ10 configurations.5An explicit calculation with τ1=τ2=0gives α=4τμ2,β=μ1−μ0.(85) 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 m2C=m2C−m2˜ C(α,β, 0)(86) produces numerically stable results with αand βgiven in (85) when the same sequence of σ0and σ1values are used in both Cand ˜ C. It turns out that the configurations of type (b) of Care also cured by the subtraction in (86). Thus, we are only left with the discussion of the case (b) for D0.Intheτ→0regionit 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, D0 iπ2=1 6m41+τ+χ 10 +τ2+χ2 70 +τχ 140 +τ3+χ3 420 +τχ(τ +χ) 1260 +Oτ4.(87) 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. Results Here we present the numerical outcome of a MC based on (71), (72), (86) and (87). The results for C(P2,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 ΔR,I:= 1+MC −AN AN R,I ,(88) where Rand Irefer to the real and imaginary parts, respectively. With the given MC statistics, the analytic result is 5Additionally, 0can be used to control when such a local subtraction has to be performed. 123
55 Page 16 of 22 Eur. Phys. J. C (2022) 82:55 Fig. 16 The one-loop self-energy diagram of (116) This gives, by definition of FDR integration, BFDR := lim μ→0ˆd4qm2 0 ¯q2¯ D0¯ D1+M2 1(q) ¯q4¯ D1μ=μR ,(119) where μRis the finite renormalization scale. It is convenient to take the limit μ→0 directly at the integrand level and substitute μwith μRonly in the logarithms. This is achieved by rewriting BFDR =ˆd4q1 D0D1−1 (q2−μ2 R+i)2,(120) where Di:= ¯ Di+μ2. By doing so, μ2is 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. 3and 4can be adapted to deal with (120). More complex multi-loop UV divergent configurations can be treated likewise. 5.1 Integrating over the loop energy component We take m1=m0=mfor simplicity. A rescaling as in (47) produces BFDR =4πˆ∞ 0 ρ2dρIt,(121) where It:= ˆ+∞ −∞ dt 1 (t2−R2 0+i)((t−√τ)2−R2 0+i) −1 (t2−R2 ν+i)2,(122) with R2 0:= ρ2+1,R2 ν:= ρ2+ν, ν := μ2 R/m2.(123) The Cauchy integral theorem allows one to compute It=iπ 21 R0−iR0−√τ/2−iR0+√τ/2−i −1 (Rν−i)3.(124) Inserting this in (121) and using r=R0−√τ/2asanew integration variable gives BFDR =2iπ2− ˆdr F(r) r−i,(125) where F(r)=θr−1+√τ/2r+√τ/22−1 ×1 r+√τ−rr+√τ/2 r+√τ/22−1+ν3/2. (126) If √τ>2 the pole at r=imigrates 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 AFDR =ˆ[d4q]1 ¯ D0 .(127) 5.2 Gluing substructures The gluing approach of Sect. 4can be easily extended to (116). Inserting (63)in(120)gives BFDR(ν) =π 2τ− ˆ1 j=0dσj σj+iλ1 2θ(λ) J(ν), (128) with J(ν) =1−σ0σ1 (σ0+μ0−ν+i)2(129) 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, BFDR(−iν) is better suited than BFDR(ν) to be evaluated numerically if ν∈ . Besides, the connection between the two can be derived by differentiating (120), ∂BFDR ∂μ2 R=−2ˆd4q (q2−μ2 R+i)3=iπ2 μ2 R ,(130) which gives BFDR(ν) =BFDR(−iν)−π3 2−iπ2ln ν ν.(131) BFDR(−iν)still suffers from numerical inaccuracies of type (76) (a). To cure this, we locally subtract from it an approximant, ˜ BFDR(−iν)=0, constructed in such a way that, after changing variables as in (79), all cuts and poles lie in the lower σ0complex half-plane. This is obtained by replacing 123
Eur. Phys. J. C (2022) 82:55 Page 17 of 22 55 Table 9 The integral of (116) with m=m0=m1=2μRas a function of √τ. The numerical estimates are computed by sampling (125)with10 8MC shots and the analytic results are obtained with OneLOop.MC errors between parentheses √τNumerical result Analytic result 0.2 2(2)×10−6−i1.3614(3)×1010−i1.3616×101 0.4 2(2)×10−6−i1.3411(3)×1010−i1.3415×101 0.6 2(2)×10−6−i1.3065(3)×1010−i1.3068×101 1.5 2(2)×10−6−i8.705(1) 0 −i8.7063 1.9 4(4)×10−7−i2.0737(3) 0 −i2.0739 2.1 −9.455(1) +i4.162(1) −9.4541 +i4.1616 4−2.6854(3)×101−i1.6453(5)×101−2.6852×101−i1.6456×101 10 −3.0377(5)×101−i3.828(2)×101−3.0380×101−i3.8280×101 50 −3.0990(6)×101−i7.112(5)×101−3.0981×101−i7.1094×101 100 −3.1009(6)×101−i8.473(7)×101−3.1000×101−i8.4825×101 Table 10 The integral of (131)withν=1/4andμ0=μ1=1. The estimatesareobtainedbysamplingwith10 8MC shots (132) evaluated at ν=1. MC errors between parentheses √τBFDR(−iν)−π3 2−iπ2ln ν ν 0.2 −4(4)×10−3−i1.357(3)×101 0.4 −3(3)×10−2−i1.342(2)×101 0.6 −1(1)×10−2−i1.308(2)×101 1.5 −1(1)×10−2−i8.69(2) 1.9 2(2)×10−3−i2.04(2) 2.1 −9.45(1) +i4.19(2) 4−2.688(2)×101−i1.655(3)×101 10 −3.039(4)×101−i3.825(5)×101 50 −3.100(8)×101−i7.12(1)×101 100 −3.09(1)×101−i8.49(2)×101 in (128)λ1 2θ(λ) →√λ2−4τiθ(λ+4τσ0). In summary, we rewrite BFDR(−iν)=π 2τ− ˆ1 j=0dσj σj+i ×λ1 2θ(λ)−λ2−4τiθ(λ+4τσ0)J(−iν). (132) In Table 10 we present our estimates for BFDR(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 108MC points, BFDR(1/4)=0.02(2)−i27.20(3), (133) to be compared to the analytic value −i27.220. 6 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 CIR := lim μ→0ˆd4q1 D0D1D2 , D0:= q2−μ2+i, D1:= (q+p1)2−μ2+i, D2:= (q−p2)2−μ2+i, (134) that appears in a P→p1+p2decay with p2 i=0 and s:= P2. 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ϕ3theory 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 xdirection. The cuts contributing to ϕ∗→ϕϕ(ϕ) are listed in Fig. 18 where, to make contact with (134), Da+De=4π2g4π 2(iC IR), Dc+Dg=4π2g4π 2(iC IR)∗.(135) The diagrams are organized in pairs sharing collinear singularities. For instance, in Dathe energy component of propagator 1 is the sum of those of propagators 2 and 3. Thus, propagators 1 and 3 never pinch in the q0 1complex plane, and propagator 3 can only become collinear to 4. Likewise, in Db 123
55 Page 18 of 22 Eur. Phys. J. C (2022) 82:55 Fig. 17 The Feynman rules of the gϕ3theory. A special notation is used for propagators with positive and negative values of p0and px, which, in our convention, coincides with the direction of the back-toback final-state particles in the Prest-frame. The complex conjugate of such rules is used in the r.h.s. of diagrams cut by a dashed line Fig. 18 The twoand three-particle cuts contributing to ϕ∗→ϕϕ(ϕ) in the Prest 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 xonly 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, 12 Note that adding μ2also to 1 and/or 2 does not change the asymptotic μ→0 limit of the result. This is used, for instance, in (135). Da+Dbis free of collinear divergences, and the same happens for Dc+Dd. In addition, Da+Db+Dc+Ddis also free of infrared singularities. A similar reasoning applies to De+Df+Dg+Dh, but with an opposite sign of the energy component of q1. The previous analysis shows that the threeparticle cuts Dband Ddcan be used as local countertems for Da+Dc.13 This requires common reference frames. One can employ two different routings for Da+Dband Dc+Dc.However, they must coincide in the limit q1→0 to guarantee the cancellation of the soft behaviour of Da+Db+Dc+Dd.In particular, when computing Da,bwe assign a momentum q2 to propagator 4, from left to right, and choose ω1:= q1/√s=(t1,ρ 1,0,0), ω2:= q2/√s=(t2,ρ 2cθ,ρ 2sθsφ,sθcφ). (136) On the other hand, we calculate Dc,dwith q2assigned to propagator 5 and ω1=(t1,−ρ1,0,0), ω2=(t2,ρ 2cθ,ρ 2sθsφ,sθcφ). (137) The result of the computation is reported in Appendix C in terms of integrals over Ri:= ρ2 i+η, with η:= μ2/s.(138) It is convenient to further split Da,c=Ds a,c+Du a,c, where the superscripts s,urefer to the subtracted and unsubtracted regions, which correspond to the integration intervals √η< R1<1/2 and 1/2<R1<∞, respectively. In fact, Db,d contribute in the subtracted region only, and Du a+Du cis free of IR singularities, s g4(Du a+Du c)=−2π5ln22+2Li2(−1/2) =254.838137 ···.(139) An analytic calculation [19] shows that j=a,b,c,dDj=0. Hence, one must have K:= s g4Ds a+Ds c+Db+Dd=−254.838137 ···. (140) In Table 11 we display our numerical estimate of Kbased 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]. 13 An analogous procedure holds for the last four cuts of Fig. 18. 123
Eur. Phys. J. C (2022) 82:55 Page 19 of 22 55 Table 11 The combination of cut-diagrams defined in (140)as a function of ηin (138). Numbers obtained with 1010 MC shots ηK 10−10 −254.81(1) 10−11 −254.83(1) 10−12 −254.84(1) Table 12 Time to generate 106MC shots on a single 2.2 GHz processor Type of integral Location Time [s] One-loop triangle Last row of Table 20.25 Two-loop self-energy Eq. (96)4.7 Two-loop vertex Last row of Table 616 Planar double box Last row of Table 871 Three-loop planar box Eq.(110)22 Two-loop pentabox Eq.(115)17 UV one-loop bubble Last row of Table 90.17 IR case Table 11 0.13 7 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 iprescription 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 109Monte Carlo shots with ∼10−7(in terms of some relevant mass scale) can yield relative precision of the order of 10−4for one-loop diagrams and 10−3for twoand 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 106MC 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. Acknowledgements The work of RP is supported by the SRA grant PID2019-106087GB-C21 (10.13039/501100011033), by the Junta de Andalucía grants A-FQM-467-UGR18 and P18-FR-4314 (FEDER), and by the COST Action CA16201 PARTICLEFACE. The work of BW was partially supported by STFC HEP consolidated grants ST/P000681/1 and ST/T000694/1. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors’ comment: This is a theoretical paper with no external associated data.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP3. Appendix A: Taylor expansions The Taylor expansions for integrals like (71) and (72) can be obtained using the general result − ˆd4qqμ1...qμkqν1...qνk (q2−m2+i)n =Cn,k iπ2 m2n−2k−4g{μ1ν1...gμkνk}(A.1) 123
55 Page 20 of 22 Eur. Phys. J. C (2022) 82:55 where Cn,k=(−1)n(−4)k(n−k−3)! k!(n−1)!,(A.2) which can be established by induction. For (71), we first make a shift of variable q→q+p2. Then for τ1=τ2=0, μ1=μ2=1wehave C=∞ k,l=0− ˆd4q(2q·p1)k(−2q·p2)l (q2−m2+i)k+l+3.(A.3) Applying (A.1), and noting that terms with l= kvanish as they would involve τ1or τ2, C iπ2=∞ k=0 (−4)kC2k+3,k m2k+2pμ1 1...pμk 1pν1 2...pνk 2 ×g{μ1ν1...gμkνk} =−∞ k=0 pμ1 1...pμk 1pν1 2...pνk 2 (2k+2)!m2k+2g{μ1ν1...gμkνk} =−1 m2 ∞ k=0 (k!)2 (2k+2)!τk,(A.4) where τ=2p1·p2/m2.Forτ1,τ 2= 0, we have instead of (A.3) C=∞ k,l=0− ˆd4q(2q·p1−p2 1)k(−2q·p2−p2 2)l (q2−m2+i)k+l+3 =−iπ2 2m21+1 12(τ +τ1+τ2)+1 90(τ2+τ2 1+τ2 2 +ττ1+ττ2+τ1τ2)+... .(A.5) This expansion can be extended to general masses by the substitutions τ1→τ1+μ2−μ1, τ2→τ2+μ2−1.(A.6) For the expansion of (72)wehave D0 iπ2= j,k,l,n l!(−m2τ) l−n n!(l−n)! ×− ˆd4q(−2p1.q)j(2p2.q)k[2(p3−p1).q]n (q2−m2+i)j+k+l+4.(A.7) Note that the number of factors of qμmust be even, j+k+ n=2K. Then D0 iπ2= j,k,l,n (−1)j+l+n4Kl!τl−n n!(l−n)! CN,K m2K+4 ×pμ1 1...pμj 1pμj+1 2...pμj+k 2 ×(p3−p1)μj+k+1...(p3−p1)μ2K ×g{μ1μ2...gμ2K−1μ2K}(A.8) where N=j+k+l+4. Substituting (A.2), we obtain (87). Appendix B: The one-point FDR integral We compute the one-loop integral AFDR =ˆ[d4q]1 ¯ D0 ,(B.9) with ¯ D0given in Eq. (117). Extracting the vacuum produces the expansion 1 ¯ D0=1 ¯q2+m2 0 ¯q4+m4 0 ¯q4¯ D0 ,(B.10) hence AFDR =lim μ→0ˆd4qm4 0 ¯q4¯ D0μ=μR .(B.11) By taking the μ→0 limit at the integrand level and replacing μwith μRin the logarithmic divergent vacuum one obtains AFDR = ˆd4q1 q2−m2 0+i−1 q2+i−m2 0 (q2−μ2 R+i)2. (B.12) Choosing now m=m0gives AFDR/m2=4πˆ∞ 0 ρ2dρIt,(B.13) where It:= ˆ+∞ −∞ dt 1 t2−R2 0+i−1 t2−ρ2+i −1 (t2−R2 ν+i)2,(B.14) with R0and Rνdefined in (123). One computes It=iπ1 ρ−1 R0−1 2R3 ν,(B.15) which gives AFDR/m2=4iπ2− ˆdρ ρθ(ρ)F(ρ), (B.16) where F(ρ) := 1 1+1 ρ21+1+1 ρ2−1 21+ν ρ23 2 . 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 AFDR/m2=iπ2(1+ln ν). (B.17) 123
Eur. Phys. J. C (2022) 82:55 Page 21 of 22 55 Table 13 The integral of (B.16) as a function of νcompared to the analytic result of (B.17). The numerical estimates are obtained with 108shots. MC errors between parentheses νNumerical result Analytic result 0.1 −i1.2855(2)×101−i1.2856×101 0.5 i3.0281(4) i3.0285 1i9.869(1) i9.8696 2i1.6709(2)×101i1.6711×101 10 i3.2596(4)×101i3.2595×101 Appendix C: The IR integrals Here we obtain onefold integral representations for the cutdiagrams Da,b,c,dof Sect. 6and prove (139). The diagram Db: Choosing the momenta as in (136)gives Db g4=8π3 sˆd4ω1d4ω2 δ+(σ2)δ+(σ3−η)δ+(σ4−η) (σ1−η+i)(σ5−i) ×θ(ρ1+ρ2cθ), (C.18) with σ1=(1−t2)2−ρ2 2, σ2=(1−t1−t2)2−ρ2 1−ρ2 2−2ρ1ρ2cθ, σ3=t2 1−ρ2 1,σ 4=t2 2−ρ2 2, σ5=σ2−1+2(t1+t2). (C.19) Note that a harmless μ2has been added to propagator 1 and that the Heaviside function forces propagator 5 to have a positive component of the momentum along x.Usingthe three Dirac delta functions one arrives at Db g4=−8π5 sˆ∞ √η dR 1ˆ∞ √η dR 2 1 1−2R2 1 1−2R+ ×θ(1−R+)θ(1−|cθ|)θ(ρ1+ρ2cθ), (C.20) with R1,2in (138), R+:= R1+R2and cθ=1 2ρ1ρ21+2(η +R1R2−R+).(C.21) Integrating analytically over R2produces logarithms with boundaries determined by the three Heaviside functions. The result reads Db g4=−2π5 sˆ1 2 √η dR 1 R1 ln ⎛ ⎜ ⎝η 1+1−2R1 R1+R2 1−η R1+R2 1−η−η⎞ ⎟ ⎠. (C.22) The diagram Da: We split Dainto two components, Da=D+ a+D− a, with positive and negative values of q1along x. Choosing the momenta as in (136) produces D+ a g4=4iπ2 sˆd4ω1d4ω2 j=1,3,41 σj−η+i ×δ+(σ2)δ+(σ5)θ(ρ1+ρ2cθ), (C.23) with the same σ1÷5of (C.19). Using the two delta functions gives D+ a g4=8iπ4 sˆ∞ 0 dρ1ρ1ˆ∞ 0 dρ2ρ2I1 ×θ(1−|cθ|)θ(ρ2 1−ρ2 2+1/4), (C.24) where cθ=1 2ρ1ρ21 4−ρ2 1−ρ2 2(C.25) and I1:= ˆ∞ 0 dt11 t2 1−R2 1+i ×1 (1/2+t1)2−R2 2+i 1 (1/2−t1)2−R2 2+i. (C.26) D− ais obtained from (C.24) by replacing θ(ρ2 1−ρ2 2+1/4)→θ(−ρ2 1+ρ2 2−1/4). Hence Da g4=8iπ4 sˆ∞ 0 dρ1ρ1ˆ∞ 0 dρ2ρ2I1θ(1−|cθ|). (C.27) Computing I1with the Cauchy integral theorem gives Da g4=−8π5 sˆ∞ √η dR 1ˆ∞ √η dR 2θ(1−|cθ|) ×1 1+2R2 1 1+2R+−1 1−2R2+i 1 1−2R++i. (C.28) Note the appearance of the same denominator structures of (C.20). An integration over R2produces Da g4=−2π5 sˆ∞ √η dR 1 R1ln 1+2R2 1+2R+ −ln 1−2R2+i 1−2R++iR+ 2 R− 2 ,(C.29) where R± 2:= $1 2±R2 1−η2 +η. (C.30) The diagram Dc: It is the complex conjugate of (C.29). 123
55 Page 22 of 22 Eur. Phys. J. C (2022) 82:55 The diagram Dd: Inserting a harmless μ2in propagator 5 and choosing the momenta as in (137)gives Dd g4=8π3 sˆd4ω1d4ω2 δ+(σ2−η)δ+(σ3−η)δ+(σ4) (σ1+i)(σ5−η−i) ×θ(cθ), (C.31) where σ1=σ4+1−2(t2−t1), σ2=(1−t2)2−ρ2 2, σ3=t2 1−ρ2 1,σ 4=(t2−t1)2−ρ2 1−ρ2 2−2ρ1ρ2cθ, σ5=t2 2−ρ2 2.(C.32) Using the delta functions produces Dd g4=−8π5 sˆ∞ √η dR 1ˆ∞ √η dR 2 1 1−2R2 1 1−2R+ ×θ(1−R+)θ(cθ)θ(1−cθ), (C.33) with cθas in (C.21). Integrating over R2gives Dd g4=−2π5 sˆR+ 1 √η dR 1 R1 L(R1), (C.34) where L(R1):= ln 1 1−2R1+1 R1+R2 1−η ×η(R1−2η) (R1+R2 1−η−η)(R1(1−2R1)+2η), (C.35) and R+ 1:= 1 2 1+4η2 1+η−4η2(1−η) .(C.36) The combination Du a+Du c: This is obtained from (C.29) and its complex conjugate by replacing the integration region √η<R1<∞by 1/2< R1<∞. Upon doing so, ηcan be set to 0 and R± 2=R1± 1/2. Introducing x=2R1gives Du a+Du c g4=−4π5 sˆ∞ 1 dx xln x+2 x+1+ln x−2 x−1, from which (139) follows. References 1. G. Heinrich, Phys. Rept. 922, 1 (2021). https://doi.org/10.1016/j. physrep.2021.03.006 2. D.E. Soper, Phys. Rev. D 62, 014009 (2000). https://doi.org/10. 1103/PhysRevD.62.014009 3. T. Binoth, G. Heinrich, Nucl. Phys. B 585, 741 (2000). https://doi. org/10.1016/S0550-3213(00)00429-6 4. T. Binoth, J.P. Guillet, G. Heinrich, E. Pilon, C. Schubert, JHEP 10, 015 (2005). https://doi.org/10.1088/1126-6708/2005/10/015 5. Z. Nagy, D.E. Soper, Phys. Rev. D 74, 093006 (2006). https://doi. org/10.1103/PhysRevD.74.093006 6. E. de Doncker, Y. Shimizu, J. Fujimoto, F. Yuasa, Comput. Phys. Commun. 159, 145 (2004). https://doi.org/10.1016/j.cpc.2004.01. 004 7. F. Yuasa, E. de Doncker, N. Hamaguchi, T. Ishikawa, K. Kato, Y. Kurihara, J. Fujimoto, Y. Shimizu, Comput. Phys. Commun. 183, 2136 (2012). https://doi.org/10.1016/j.cpc.2012.05.018 8. E. de Doncker, F. Yuasa, K. Kato, T. Ishikawa, J. Kapenga, O. Olagbemi, Comput. Phys. Commun. 224, 164 (2018). https://doi. org/10.1016/j.cpc.2017.11.001 9. J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, M. Spira, J. Streicher, JHEP 04, 181 (2020). https://doi.org/10.1007/ JHEP04(2020)181 10. A. Ghinculov, Phys. Lett. B 385, 279 (1996). https://doi.org/10. 1016/0370-2693(96)00871-4 11. Guillet, J. Ph. and Pilon, E. and Shimizu, Y. and Zidi, M. S., PTEP 2020(4), 043B01 (2020). https://doi.org/10.1093/ptep/ptaa020 12. Bauberger, Stefan and Freitas, Ayres and Wiegand, Daniel, JHEP 01, 024 (2020). https://doi.org/10.1007/JHEP01(2020)024 13. R. Pittau, JHEP 11, 151 (2012). https://doi.org/10.1007/ JHEP11(2012)151 14. R. Kleiss, R. Pittau, Comput. Phys. Commun. 83, 141 (1994). https://doi.org/10.1016/0010-4655(94)90043-4 15. D. Kermanschah, e-Print: arXiv:2110.06869 [hep-ph] (2021) 16. A. van Hameren, Comput. Phys. Commun. 182, 2427 (2011). https://doi.org/10.1016/j.cpc.2011.06.011 17. D.J. Broadhurst, Z. Phys. C 47, 115 (1990). https://doi.org/10. 1007/BF01551921 18. G. Ossola, C.G. Papadopoulos, R. Pittau, JHEP 03, 042 (2008). https://doi.org/10.1088/1126-6708/2008/03/042 19. R. Pittau, Eur. Phys. J. C 74(1), 2686 (2014). https://doi.org/10. 1140/epjc/s10052-013-2686-1 20. C. Gnendiger et al., Eur. Phys. J. C 77(7), 471 (2017). https://doi. org/10.1140/epjc/s10052-017-5023-2 123