Abstract
An improved method is presented for the numerical evaluation of multiloop integrals in dimensional regularization. The technique is based on MellinBarnes representations, which have been used earlier to develop algorithms for the extraction of ultraviolet and infrared divergencies. The coefficients of these singularities and the nonsingular part can be integrated numerically. However, the numerical integration often does not converge for diagrams with massive propagators and physical branch cuts. In this work, several steps are proposed which substantially improve the behavior of the numerical integrals. The efficacy of the method is demonstrated by calculating several twoloop examples, some of which have not been known before.
On the Numerical Evaluation of Loop Integrals
[.5ex] With MellinBarnes Representations
[3.5em] Ayres Freitas, YiCheng Huang
Department of Physics & Astronomy, University of Pittsburgh,
3941 O’Hara St, Pittsburgh, PA 15260, USA
1 Introduction
Perturbation theory is one of the most important tools for the calculation of particle physics observables. In many cases, radiative corrections are important to match the experimental precision. General and efficient techniques for the computation of oneloop corrections have been developed over many years, resulting in automated frameworks for the calculation of processes with up to six external legs [1]. Most methods for oneloop corrections reduce the loop integrals to a welldefined set of master integrals by PassarinoVeltman and Melrose reduction [2], integrationbyparts and Lorentzinvariance relations [3], or onshell cut techniques [4]. See Ref. [5] and references therein for a more detailed overview.
In some cases, however, oneloop corrections are not sufficient and two or even threeloop corrections becomes relevant. Typical examples are electroweak precision observables [6, 7] and production of gauge bosons [8], top quarks [9, 10] or Higgs bosons [11] at the Large Hadron Collider (LHC). Many highly complex calculations have been performed for this purpose by several authors (see references in previous sentence), involving specialized techniques tailored to each particular process. For twoloop amplitudes with a small number of different mass and momentum scales one can follow the traditional approach of reduction to a small number of master integrals, which are then calculated analytically. However, this strategy does not work for applications with a large number of different scales, due to the complexity of the expressions and the fact that in general the master integrals cannot be determined analytically beyond the oneloop level.
Fully numerical techniques for the evaluation of two and higherloop integrals face two main challenges: extraction of ultraviolet (UV) and infrared and collinear (IR) singularities, as well as stability and efficiency of the numerical integrations. Methods based on dispersion relations [12, 7] and integrations over Feynman parameters [13] result in lowdimensional numerical integrals with good convergence properties, but they do not offer a general algorithm for the isolation of divergencies. Two powerful approaches for a straightforward extraction of poles in dimensional regularization are sector decomposition [14] and MellinBarnes representations [16, 17, 15]. However, these methods lead to complicated multidimensional numerical integrals, which have convergence problems for diagrams with physical branch cuts. For sector decomposition, the instability related to physical branch cuts can be circumvented by the integration contour deformation proposed in Ref. [18].
This article instead focuses on the numerical integration of MellinBarnes (MB) representations. Compared to sector decomposition, MB representations of Feynman integrals have the advantage that they typically lead to substantially shorter expressions in the integrand. However, MB representation will only become a viable general tool for numerical twoloop integrations if the behavior of the numerical integrals can be substantially improved. In this paper, several techniques for this purpose are presented, which can also be used in combination to optimize the convergence properties of the integrals. In particular, complex variable transformations and relations that reduce the number of integration dimensions prove to be very useful.
Section 2 starts out by reviewing the employment of MB representations for loop integrals and explaining the origin of the numerical instabilities in the MB integrals. In section 3 the new methods to deal with these instabilities are described in detail. Numerical results for several examples of twoloop integrals are presented in section 4. Where applicable, they are compared to existing results in the literature, which have been obtained by other methods. In addition, results are shown for some integrals appearing in the calculation of twoloop QCD corrections to which have not been calculated before. Finally, the main findings are summarized in section 5.
2 MellinBarnes representations
The following discussion of the construction of MB representations and the isolation of UV and IR singularities is largely based on Ref. [16]. Let us start with a oneloop point integral
(1) 
and its Feynman parametrization
(2)  
where . The integrand denominator can be transformed into a MB representation by
(3)  
where the following two conditions must be met:

The integration contours for are straight lines parallel to the imaginary axis chosen such that all arguments of the gamma functions have positive real parts.

The are complex numbers with for any pair .
Satisfying the second condition potentially requires some adjustments of the Feynman parametrization. For example, the scalar diagram in Fig. 1 is represented by the Feynman integral
(4) 
where the terms with violate condition 2 above. More generally, problems with the application of the MB representation arise if one mass or momentum parameter appears both with positive and negative signs in the denominator of the Feynman integral. This complication can be avoided by using the relation . In the case of eq. (4), this relation allows to rewrite it as
(5) 
so that all terms with now have positive signs. The term does not lead to any difficulties since one can imagine to temporarily shift slightly off the real axis into the complex plane so that is satisfied. Thus the MellinBarnes formula can be safely applied to eq. (5).
After the introduction of the MB representation, the integration over the Feynman parameters can be carried out straightforwardly, using
(6) 
assuming that the exponents satisfy
(7) 
A MB representation for a twoloop integral can be obtained by first deriving a MB representation for one subloop, as described above, and inserting the result into the second loop, which then can be transformed into a MB representation with the same procedure. In general, the MB representation for the first subloop will contain invariants depending on the second loop momentum, which can be interpreted as new propagators for the second integration.
As an example, consider the “sunset” diagram in Fig. 2. The MB representation for the upper subloop is given by
(8) 
where and is the internal momentum of the second (lower) loop. For the second integration one then obtains
(9) 
so that the gives rise to a new propagator. After constructing the MB representation for the second loop, the final result is
(10) 
The MB integral can then be be performed numerically.
As mentioned above, the MellinBarnes formula and the Feynman parameter integration require the arguments of all gamma functions to have a positive real part. In general, this requires to differ from zero by a finite amount. When taking the limit , the poles of some gamma functions can cross some of the integration contours, so that one needs to correct for this by picking up the residue of the integrand for crossed each pole. As a result, for one ends up with the original integral plus a sum of terms with lower integration dimension stemming from these residues. An algorithm for this procedure is described in Refs. [16, 17].
Since for appropriate choices of the contours each integral in this sum is regular, the singularities must be contained in the coefficients of these integrals. Therefore this method allows to extract all UV and IR divergencies in a systematic way.
3 Improving numerical convergence
The contour integrations in the MB representation can be parametrized as
(11) 
where the are real constants. The integrand contains powers of the kinematic and mass invariants as well as gamma functions.
Let us return to the example of the sunset diagram, eq. (10). For Euclidean momentum, , the dependence of the invariants on the integration variables is given by simple oscillating exponentials, , while the gamma functions rapidly tend to zero when their arguments take values away from the real axis. Therefore, in this case, the integrals are well behaved for numerical integration, and a moderate finite integration interval is sufficient to reach a high precision for the result.
However, for physical momentum, , an exponentially growing term appears in the integrand:
(12) 
Thus the integrand is an oscillating function with slowly decreasing amplitude for . While the integral is formally convergent, numerical integration routines have difficulties dealing with this situation, since one needs to integrate to very large values of , over many oscillation periods, to reach acceptable precision. See Fig. 3 for illustration.
The problem of slowly converging oscillating tails of the integrand is generic for any MB representation of a loop integral with physical (Minkowski) momenta. Special techniques for numerical integration of oscillating integrands like Filon’s method [19] fail since the oscillation period varies substantially within the integration range. Instead, as shown in the following, suitable variable transformations can substantially improve the convergence properties.
Next, one can exploit the freedom to deform the integration contours in the complex plane, as long as no contour crosses a pole of the integrand. The latter condition is rather restrictive, but if all contours are rotated in the complex plane by the same angle it is guaranteed that no pole crossing occurs:
(13) 
In eq. (12) this leads to
(14) 
so that the problematic exponential can, in principle, be compensated by the a suitable choice of . For the example shown in Fig. 3 (right), it is evident that with this contour deformation the integrand drops off very fast away from the origin, and thus is much better suited for numerical integration.
For a multidimensional MB integral, the choice of is more delicate since this one parameter affects all integration variables. It helps to replace the variables by hyperspherical coordinates, angles and one radial coordinate:
(15)  
In this representation, only is affected by the modification (13) of the integration contours, so that one can in principle choose such that all problematic exponentials of the form (12) are cancelled, as a function of the .
After applying the variable transformations, the multidimensional MB integral is sufficiently wellbehaved that it can be evaluated with standard MonteCarlo integration algorithms like Vegas [20]^{*}^{*}*For relatively simple integrals with four or less integration dimensions, deterministic integration methods like Gaussian quadrature prove more efficient.. In particular, the MonteCarlo error decreases systematically with the number of integration points, which is an important requirement for the convergence of the numerical integration routine.
However, MB integrals with high dimensionality can still require a large amount of computing time to reach a result with acceptable precision. For twoloop diagrams one can obtain MB integrals with more than 10 dimensions, depending on the topology and the number of independent masses and momenta involved. Fortunately, it turns out that in most cases some of the integrations can be performed analytically. For this purpose, the convolution theorem for Mellin transforms proves to be useful. The Mellin transform of a function is defined as
(16) 
For the Mellin transform of a product of two functions, the convolution theorem then states
(17) 
where the integration is parallel to the imaginary axis and the integrand in eq. (17) must be analytic for Re in some region around .
With suitable choices of and one can derive formulas for MB integrals of two or four gamma functions:
(18)  
(19)  
(20)  
where in each case the contours are assumed to be straight lines parallel to the imaginary axis satisfying the conditions in square brackets, and is the Gauss hypergeometric function. These formulas should be applied before the variable transformations. The maximal reduction of the integration dimension is achieved if the formulas are used both before and after the extraction of the poles, since in each case one obtains different combinations of terms that match the lefthand sides of eqs. (18)–(20). The reduction procedure has been implemented in a computer algebra program based on Mathematica, which automatically looks for all possible applications of these equations.
4 Examples
In this section numerical results for some examples of scalar twoloop integrals are presented. The generation of the MB representation, the extraction of the poles, and the application of the reduction formulas, as described in the previous section, has been implemented in Mathematica. For the numerical integration a C code has been developed, which uses Gaussian quadrature for low dimensionality and the Vegas MonteCarlo routine for high dimensionality.
In order to test the method and its computer implementation, results will be shown first for several known twoloop integrals, which have been calculated earlier with other methods. A few examples that cannot be computed with any previously published method will be presented at the end of this section.
It should be stressed that all examples in this section cannot be evaluated numerically with MB representations by using the standard representation eq. (11), as implemented in the mb package of Ref. [17], see also Ref. [16]. The reason is that for physical momenta, , the integrands have long oscillating tails that converge too slowly. Therefore, for all cases shown in the following, the variable transformations and contour deformations introduced in the previous sections have been used.
Selfenergy topologies:
The scalar sunset diagram in Fig. 2 leads to a threedimensional MB integral, as shown in section 2, together with one and twodimensional integrals originating from pole crossings when the limit is taken. Expanding in powers of , it turns out that the UV poles are given by terms without integrals:
(21) 
Numerical results for the finite part are shown in table 1 for two choices of parameters, comparing with values obtained through the dispersion relation method of Ref. [12].
(1, 1, 4, 5)  (16, 2, 2, 1)  

Result from MB repr.  18.19558(2)  
(with integration error)  
Result from disp. rel. 
The scalar fivepropagator selfenergy in Fig. 4 can be cast into a tendimensional MB representation, which can be reduced to integrals with dimension up to six with the help of the reduction formulas presented in the previous section. Note that this integral is UVfinite, so that no terms from pole crossings occur.
Owing to the relative large dimensionality of the integration, the numerical integration is much slower than in the previous example. With about one hour evaluation time on a single core of a Pentium IV computer with 2.4 GHz a precision of 0.1% to 1% is obtained, depending on the values of the masses and the momentum, see table 2. They compare well with results based on the dispersion relation method of Ref. [12].
(, 1, 1, 1, 1, 1)  (1, 3.7, 3.7, 1, 3.7, 3.7)  (20, 1, 1, 1, 1, 1)  

Result from MB repr.  
(with integration error)  
Result from disp. rel. 
Vertex topologies:
Figure 5 shows a few examples for twoloop vertex contributions with three or four different scales. The three scalar diagrams lead to MB integrals with 9, 10, and 11 dimensions, respectively, which can be reduced to five dimensions each using the reduction formulas in eqs. (18) ff. In table 3, numerical results are compared with numbers obtained with the dispersion relation method of Refs. [12, 7].
Diagram  Fig. 5 (a)  Fig. 5 (b)  Fig. 5 (c) 

Result from MB repr.  
(with integration error)  
Result from disp. rel. 
Box topologies:
This subsection will discuss some twoloop box integrals corresponding to typical topologies for nexttonexttoleading order QCD corrections to , shown in Fig. 6. In the diagrams massless (thin) lines correspond to light quark or gluon propagators, while massive (thick) lines indicate top quark propagators.
The master integrals corresponding to (a) and (b) in the figure lead to MB integrals with dimension up to two and three, respectively, after expansion in . Therefore they can be evaluated very efficiently with deterministic integration methods, like Gaussian quadrature. These master integrals were computed previously in Ref. [10] with the method of differential equations. Comparing results between the two methods^{†}^{†}†We are grateful to the authors of Ref. [10] for kindly providing numerical results of their calculation., agreement to more than nine significant digits has been obtained^{‡}^{‡}‡Higher precision could be achieved with more computing time.. This comparison was carried out both for an unphysical parameter point with Euclidean momentum , as well as for a physical parameter point with . For these cases, the evaluation of the MB integrals takes only several seconds to minutes on a singlecore Pentium IV computer.
To the best of our knowledge, the twoloop integral corresponding to Fig. 6 (c) and (d) have not been computed before. Fig. 6 (c) leads to a MB representation with a 13dimensional integral, which can be reduced to six dimensions with the help of the reduction formulas. Numerical results for a few parameter points are shown in table 4. Since the structure of the integral is rather complicated, the evaluation is much slower than the previous examples. On a singlecore Pentium IV computer with 2.4 GHz, the numerical evaluation of one parameter point takes 45 hours to reach the precision reported in the table.
Result from MB representation (with integration error)  

5 Summary
MellinBarnes representations are a powerful tool for the calculation of multiloop integrals. They have been shown to be very useful for systematically extracting poles of the form in dimensional regularization. The coefficients of these poles and the finite part can then be evaluated as numerical integrals. However, for diagrams with massive propagators and physical branch cuts, the integrands a priori can be highly oscillatory, so that standard numerical integration routines do not yield convergent results.
In this article, it is shown how the numerical convergence behavior can be improved substantially by transforming the integration variables, deforming the integration contours in the complex plane, and using analytical formulas to reduce the dimensionality of the integration. The usefulness of those modifications was demonstrated by calculating several examples of twoloop selfenergy, vertex and box topologies, some of which have been computed for the first time.
Several previously known techniques permit much faster evaluation of relatively simple multiloop integrals with few independent masses and momenta than the method presented here. However, MellinBarnes representations can be applied in principle to any loop integral, irrespective of its complexity. With sufficient computing time, arbitrarily precise results can be obtained.
In this work, tensor integrals with external uncontracted Lorentz indices were not discussed separately. In fact, for most applications, it is convenient to decompose the loop amplitudes into standard matrix elements or form factors, the coefficients of which only contain scalar integrals. The integrands of these scalar integrals can contain scalar products in the numerator involving on one or more of the loop momenta. Such integrals can be evaluated straightforwardly with MellinBarnes representations by expressing the scalar products as propagators with negative powers.
Acknowledgements
The authors would like to thank R. Bonciani, A. Ferroglia, and C. Studerus for useful discussions and providing numbers for comparison. This work has been supported in part by the National Science Foundation under grant no. PHY0854782.
References

[1]
A. Denner, S. Dittmaier, M. Roth and L. H. Wieders,
Phys. Lett. B 612, 223 (2005);
A. Denner, S. Dittmaier, M. Roth and L. H. Wieders, Nucl. Phys. B 724, 247 (2005);
T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and T. Reiter, Comput. Phys. Commun. 180, 2317 (2009);
R. K. Ellis, W. Giele, Z. Kunszt, K. Melnikov, G. Zanderighi, JHEP 0901, 012 (2009);
R. K. Ellis, K. Melnikov and G. Zanderighi, JHEP 0904, 077 (2009);
C. F. Berger et al., Phys. Rev. Lett. 102, 222001 (2009);
C. F. Berger et al., Phys. Rev. D 80, 074036 (2009). 
[2]
D. B. Melrose,
Nuovo Cim. 40, 181 (1965);
G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160, 151 (1979). 
[3]
F. V. Tkachov,
Phys. Lett. B 100, 65 (1981);
K. G. Chetyrkin and F. V. Tkachov, Nucl. Phys. B 192, 159 (1981);
T. Gehrmann, E. Remiddi, Nucl. Phys. B 580, 485 (2000). 
[4]
Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower,
Nucl. Phys. B 425, 217 (1994);
Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 435, 59 (1995);
Z. Bern, L. J. Dixon and D. A. Kosower, Phys. Rev. D 73, 065013 (2006);
G. Ossola, C. G. Papadopoulos and R. Pittau, Nucl. Phys. B 763, 147 (2007);
C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, JHEP 0703, 111 (2007);
R. K. Ellis, W. T. Giele and Z. Kunszt, JHEP 0803, 003 (2008);
D. Forde, Phys. Rev. D 75, 125019 (2007);
W. T. Giele, Z. Kunszt and K. Melnikov, JHEP 0804, 049 (2008);
C. F. Berger et al., Phys. Rev. D 78, 036003 (2008).  [5] Z. Bern et al. [NLO Multileg Working Group], in Proc. of the 5th Les Houches Workshop on Physics at TeV Colliders, Les Houches, France, 1129 Jun 2007, arXiv:0803.0494 [hepph].

[6]
M. Awramik, M. Czakon, A. Freitas and G. Weiglein,
Phys. Rev. D 69 (2004) 053006;
M. Awramik, M. Czakon, A. Freitas, G. Weiglein, Phys. Rev. Lett. 93, 201805 (2004);
M. Awramik, M. Czakon and A. Freitas, Phys. Lett. B 642, 563 (2006);
W. Hollik, U. Meier and S. Uccirati, Nucl. Phys. B 731, 213 (2005);
W. Hollik, U. Meier and S. Uccirati, Nucl. Phys. B 765, 154 (2007);
M. Awramik, M. Czakon, A. Freitas and B. A. Kniehl, Nucl. Phys. B 813, 174 (2009).  [7] M. Awramik, M. Czakon and A. Freitas, JHEP 0611, 048 (2006).

[8]
R. Hamberg, W. L. van Neerven and T. Matsuura,
Nucl. Phys. B 359, 343 (1991)
[Erratumibid. B 644, 403 (2002)];
C. Anastasiou, L. Dixon, K. Melnikov and F. Petriello, Phys. Rev. D 69, 094008 (2004);
K. Melnikov and F. Petriello, Phys. Rev. D 74, 114017 (2006). 
[9]
M. Czakon, A. Mitov and S. Moch,
Phys. Lett. B 651, 147 (2007);
M. Czakon, A. Mitov and S. Moch, Nucl. Phys. B 798, 210 (2008);
M. Czakon, Phys. Lett. B 664, 307 (2008);
R. Bonciani, A. Ferroglia, T. Gehrmann, D. Maître and C. Studerus, JHEP 0807, 129 (2008).  [10] R. Bonciani, A. Ferroglia, T. Gehrmann and C. Studerus, JHEP 0908, 067 (2009).

[11]
R. V. Harlander and W. B. Kilgore,
Phys. Rev. Lett. 88, 201801 (2002);
C. Anastasiou and K. Melnikov, Nucl. Phys. B 646, 220 (2002);
C. Anastasiou, K. Melnikov and F. Petriello, Phys. Rev. Lett. 93, 262002 (2004);
G. Davatz et al., JHEP 0607, 037 (2006);
S. Catani and M. Grazzini, Phys. Rev. Lett. 98, 222002 (2007);
S. Actis, G. Passarino, C. Sturm and S. Uccirati, Phys. Lett. B 670, 12 (2008). 
[12]
S. Bauberger, F. A. Berends, M. Böhm and M. Buza,
Nucl. Phys. B 434, 383 (1995);
B. A. Kniehl, Acta Phys. Polon. B 27, 3631 (1996). 
[13]
P. Cvitanovic and T. Kinoshita,
Phys. Rev. D 10, 4007 (1974);
A. Ghinculov and J. J. van der Bij, Nucl. Phys. B 436, 30 (1995);
T. Kinoshita and M. Nio, Phys. Rev. D 70, 113001 (2004). 
[14]
T. Binoth and G. Heinrich,
Nucl. Phys. B 585, 741 (2000);
T. Binoth and G. Heinrich, Nucl. Phys. B 680, 375 (2004). 
[15]
V. A. Smirnov,
Phys. Lett. B 460, 397 (1999);
J. B. Tausk, Phys. Lett. B 469, 225 (1999).  [16] C. Anastasiou and A. Daleo, JHEP 0610, 031 (2006).
 [17] M. Czakon, Comput. Phys. Commun. 175, 559 (2006).

[18]
Z. Nagy and D. E. Soper,
Phys. Rev. D 74, 093006 (2006);
C. Anastasiou, S. Beerli and A. Daleo, JHEP 0705, 071 (2007).  [19] L. N. G. Filon, Proc. Roy. Soc. Edin. 49, 38 (1928).
 [20] G. P. Lepage, “Vegas: An Adaptive Multidimensional Integration Program,”, report CLNS80/447 (1980).