Numerical Evaluation of Tensor Feynman Integrals in Euclidean Kinematics

For the investigation of higher order Feynman integrals, potentially with tensor structure, it is highly desirable to have numerical methods and automated tools for dedicated, but sufficiently 'simple' numerical approaches. We elaborate two algorithms for this purpose which may be applied in the Euclidean kinematical region and in d=4-2*epsilon dimensions. One method uses Mellin-Barnes representations for the Feynman parameter representation of multi-loop Feynman integrals with arbitrary tensor rank. Our Mathematica package AMBRE has been extended for that purpose, and together with the packages MB (M. Czakon) or MBresolve (A. V. Smirnov and V. A. Smirnov) one may perform automatically a numerical evaluation of planar tensor Feynman integrals. Alternatively, one may apply sector decomposition to planar and non-planar multi-loop epsilon-expanded Feynman integrals with arbitrary tensor rank. We automatized the preparations of Feynman integrals for an immediate application of the package sector_decomposition (C. Bogner and S. Weinzierl) so that one has to give only a proper definition of propagators and numerators. The efficiency of the two implementations, based on Mellin-Barnes representations and sector decompositions, is compared. The computational packages are publicly available.


Introduction
One goal of present calculations in particle physics is reaching higher and higher precision in perturbation theory. Since Feynman's time we have rules allowing to build automatically the necessary mathematical objects. For a long term these were just the Feynman diagrams, but recently other approaches like the unitarity based perturbative approach get rising attention. The problem remains how to evaluate the complicated integrals originating in the perturbative picture of elementary particle interactions. Physics predicts these integrals to be defined in Minkowskian space-time. The integrals are multi-dimensional complex functions with a complicated singularity structure. Further, they have to be regularized by notions like d-dimensional space-time. In recent years ambitious projects appeared where techniques are used to calculate physical processes in highly automatized ways [1,2,3].
Here we focus on the calculation of Feynman integrals in Euclidean space-time. Though they are, in general, not the ultimate physical objects, their knowledge is very useful. First, if we know them analytically, analytic continuation gives a way to transform them to the Minkowskian region, if needed so. Moreover, if we solve, somehow, Feynman integrals involved in a given physical process analytically, we can use the knowledge in the Euclidean region to check these solutions numerically. This has been done in numerous cases at the 2-loop level (e.g. massive Bhabha scattering [4,5,6,7,8,9], QCD calculations [10,11,12]), but also at higher loop levels or in more general context [13,14,15,16]). For structural studies of quantum field theory, the numerical calculation of Feynman integrals in Euclidean space-time has also been proven to be useful, e.g. for the check of some conjectures in super-Yang-Mills theories [17,18].
In the present article we describe two publicly available computational tools based on Mathematica which facilitate numerical calculations of the kind described.
The first tool is the extended version v.2.0 of the AMBRE program [19] which generates Mellin-Barnes (MB) representations for Feynman integrals. We discuss the construction of MB-integrals with numerators of arbitrary rank R for L-loop cases. We also shortly report on additional features of older versions (v.1.1 and v.1.2) which were released after publication of [20]. We explicitely work out a variety of non-trivial numerical examples. For this purpose, AMBRE is being combined with the Mathematica packages MB [21] and MBresolve [22].
The second tool is the Mathematica interface CSectors to the Ginac package sector_decomposition [23]. The package sector_decomposition uses the sector decomposition method to calculate general polynomial structures which are present in calculations of Feynman integrals. In the spirit of AMBRE, we perform in CSectors for a given Feynman integral the automatic calculation of the characteristic F and U polynomials, add some normalizing factors consisting of Gamma functions and, finally, use a general formula for multi-loop tensorial Feynman parameterisations. The result is a user-friendly interface to sector_decomposition for the specific purpose of tensor Feynman integral calculations. What remains to be done by the user of CSectors/sector_decomposition is writing in a proper way the definitions of propagators (plus numerators, if they are present). Further, optional algorithmic strategies have to be chosen which are part of the core program [23]. This is the stage of automatization reached also with AMBRE.
The article is organized as follows. In Section 2 we prepare expressions for the general multi-loop Feynman integral. Their evaluation based on Mellin-Barnes representations with AMBRE is described in Section 3. In Section 4 some details on sector decomposition and of using CSectors are discussed, and Section 5 contains numerical examples and few comparisons of the two approaches. It follows the Summary. In the appendix we list the most important Mathematica functions of AMBRE and options for CSectors.

Definitions
Comprehensive overviews of the presentations of Feynman integrals may be found e.g. in [24,25].
Here, we repeat some basic formulae in order to define our notations. The L-loop Feynman integral in d = 4 − 2ε dimensions with N internal lines with momenta q i and masses m i , and E external legs with momenta p e is defined here as follows: The numerator T R (k) is a tensor of rank R in the integration variables: and We allow for arbitrary indices ν i , the powers of propagator functions D i . Next, the momentum integrals are replaced by Feynman parameter integrals: where N ν = ∑ N i=1 ν i . The two functions U and F are characteristics of the topology of the Feynman integral. One may derive them from namely: The object A r P R−r contains the tensor structure due to its two elements: and P 0 = 1, (2.16) P µ 1 ···µ r r = P µ 1 · · · P µ r , (2.17) where we left out in the notations the indices related to the loop numbering, because they are fixed by (2.2) when the Lorentz indices are defined: The product A [µ 1 ,...,µ r r P µ r+1 ,...,µ R ] R−r is completely symmetrized in its Lorentz indices; take as an example A 2 P 2 : 20) and, more explicitely, to e.g. T (k 2 ) correspond the terms: and, with a different numerical factor (see (2.5)): For one-loop integrals, the rotation matrix M 1 in the loop momenta becomes trivial, M 1 = M 1 = 1, and so also U 1 = det(M 1 ) = 1, and F 1 (x) = −J + Q 2 . 1 The (2.5) then becomes: (2.23) In case of e.g. L = 1, R = 2, we get for the sum in (2.5): The general expressions as well as the examples agree with [19]. For one-loop tensors, equivalent expressions are also given in [26]. An important observation is that the Feynman parameter representations for arbitrary Feynman integrals would be just dependent on polynomials in the x i , i.e. be sums of monomials in the x i with integer exponents, if there were not the two types of terms U(x) A(L,N,d,R) and F(x) B(L,N,d,r) . The functions U(x) and F(x) are such polynomials, but they have non-integer exponents. While, the additional terms arising from the tensorial structure of the L-loop Feynman integrals are polynomials in the x i . 2 The function U(x) = ∑ n m n (x i ) depends only on monomials m n (x i ) and is positive semi-definite, while function F(x) = ∑ n [−s n ]m n (x i ) + U(x) ∑ N j x j m 2 j depends also on the kinematical invariants and on the masses of the problem. For Euclidean kinematics, all the [−s n ] ≥ 0 and also F(x) becomes positive semi-definite [25]. One typical example is the F-function (3.4).
The above formulae may be used for automated evaluations of specific Feynman integrals. For that purpose, one has to develop methods for their proper treatment, and two of them are worked out here.

Integrations using Mellin-Barnes representations
may be used to transform the x-integrand of (2.5) into a sequence of monomials in the x i , allowing thus to perform the x-integrations applying One remains with the problem to evaluate multi-dimensional complex Mellin-Barnes integrals. The publicly available package AMBRE prepares these Mellin-Barnes integrals. Since the first publication of AMBRE in [19], several features has been added to the package. In AMBREv.1.1, the MB-representations can be constructed only for tensor integrals where all momenta of the tensor T (k i ) are multiplied by external momenta, i.e. are part of scalar products. Since AMBREv.1.2 it is foreseen to generate MB-representations for just tensor integrals. Additionally, some new options were added, consult for details on them the webpages [20], where also appropriate examples are documented. One of the options allows to generate Feynman parameter representations without performing the x-integrations, leaving them to be performed by the preferred technique of the user.
Here we focus on AMBRE v.2.0 which generates tensor MB-integrals in a fully automatic way. Previous versions have the option Fauto which allows for manual manipulations on the F polynomial. Sometimes this is useful and helps to obtain a smaller dimensionality of MB-integrals (see e.g. [19] and the discussion for pentagons there). However, for many real processes with a large number of amplitudes, complete automatization is necessary. A second goal of the present version is a construction of MB-representations for tensorial planar n-loop tensors. Though non-planar topologies can be also obtained directly with the loop-by-loop method employed in AMBRE [27], results may come out wrong. 3 We recommend not to use AMBRE without independent checks for non-planar problems with several kinematical scales; they hopefully shall be treated more properly in the future.
In order to include tensor structures of Feynman integrals, we had to modify properly the iterative procedure of [19]. The best way to explain this might be an explicit example, see Here all propagators are assumed to be massless and the numerator is of rank R = 3: The calculation starts by working out the sub-loop integration over k 1 , which leads to the following F polynomial: According to (2.23), the generated tensor structure for the first sub-loop is an expression with five separate parts; we omit here a normalizing factor and end up with a table of five different integrals: Evidently, we will have to perform the second sub-loop integration over k 2 separately for all the above parts because they have different tensor ranks. As a rule, the rank of a given integral in the next step will include higher rank tensors than the original one. The situation after the first momentum integration can be symbolized as follows: where the expressions MB i stand for different Mellin-Barnes parts of the net integral. These parts are different due to different Feynman parameters in (3.5). As one explicite example we reproduce MB 1 : Working with more than two loops, additional iterations can produce further 'fragmentations' of the expression.
In the program, all above steps are hidden to the user. The only action is to define the input object: The output is: The dots stand for the remaining four MB-integrals. The complete expression can be found in the file MB_SE5l0m.nb at the webpage [20].

Integrations by sector decomposition
The second approach -performing sector decompositions -transforms the x-integrand into a sequence of expressions where the singularities at d = 4, as a function of ε (d = 4 − 2ε), are separated such that the arising expressions can be smoothly integrated. An appropriate algorithm for the automated computation of the ε series of multi-loop integrals has been formulated in [30,31]. In [23], a publicly available program is described which calculates the Laurent expansion of the following type of parametric integrals in ε: The a i , b i , c j and d j are integers, and the P j are polynomials in the variables x 1 , ..., x n . The program may handle a product of several polynomials, and it is not required that the polynomials are homogeneous. These features are important for applying (4.1) to our tensorial structures: it allows to calculate G L (T ) (2.5), and some examples are given in [23] for scalar integrals.
After "Generating c++ source" and "Compiling source code" in (4.4), the C++ files are running. They have the structure discussed in [23] and a user can inspect them by switching in the option "TempFileDelete->False" in (4.2). Other useful options are listed and described in the Appendix.
By default, the numerical errors of the results are also given. They are calculated by taking the combined error for all the integrals I i calculated at a given term of the ε-expansion: (4.5) As it is well known, some of integrals, especially massless ones, can be difficult to integrate numerically. Here the proper choice of a sector decomposition strategy may help [23]. In (4.4), the calculation has been done with optional strategy C. In order to control which ε-term is currently calculated, the default option TempFileDelete->False produces a log file. For the Q11 term of the example, it is: In order to test CSectors we have checked many topologies: multi-loop tadpoles and selfenergies, three-loop vertices up to rank R = 5 and with double dots on propagators (corresponding to setting index ν = 3, because a 'dot' raises the index by one), four-point functions up to rank R = 4, and some one-loop five-and six-point functions. For some higher rank tensors we have used Integration-by-Parts decompositions of the Feynman integrals using the computer algebra package IdSolver 4 . Some numerical examples of these tests can be found in the file numerical_checks.nb at the webpages [32].

Numerical results
In Section 3 it has been shown how to define Feynman integrals and how to get numerical values for them at chosen kinematical points using the MB-method. For sector decomposition, the same has been discussed in Section 4. For scalar integrals it is straightforward to use, together with AMBRE, the MB package, and to perform the numerical integrations [21]. For tensor integrals, especially with loop order L > 1, we have usually many MB-integrals for which the command MBrules of MB has to be performed in order to find proper integration paths for the MB-integrations. Sometimes, depending on the degree of divergency of the original Feynman integral, this is not possible if only an analytical continuation in ε is done. Then, an analytical continuation in one of the indices, may be successful. This has also been automatized. If MBrules does not find a valid rule for some ε = 0, then the power of the first propagator ν 1 is changed, e.g. if ν 1 = 1, then ν 1 → ν 1 = 1 + η is applied. If again MBrules cannot find a valid rule, ν 2 is treated analogously, and the procedure can be continued until MBrules is successful. 5 We introduced the auxiliary file MBnum.m to MB.m which realizes this procedure and the subsequent automatized analytic continuation, ε-expansion, and numerics. The file may be obtained from the webpages [20].
If an already prepared MB-representation repr is available, e.g. in (3.8), it is enough to use the MBnum function (see Appendix for details), e.g.: Another example is two-loop Bhabha scattering. For two-loop 4-point functions, sector decomposition needs a lot of RAM (typically up to few GB) and also of computing time. In such cases, often the numerical integrations are done faster using the MB-method. Some numerical results for the Bhabha Feynman integral with rank R = 2, corresponding to the topology shown in Figure 5.1, can be found in Table 5.1. The Mandelstam variables are s = (p 1 + p 2 ) 2 and t = (p 2 + p 4 ) 2 . The so-called reducible numerators of a tensor Feynman integral can be contracted with propagators. Apart from speeding up calculations, such contractions have been used to check the implementation of tensor structures into AMBRE and CSectors. In our example, the relation may be used to change the rank R = 2 tensor of (5.3) into three tensor integrals of rank R = 1 (in addition reducing the number of propagators by one for two of them). In this way, the numerical results given in Table 5.1 can be cross-checked.
Further, let us present a four-loop self-energy scalar diagram, Fig. 5.2. With CSectors it takes a few hours to calculate the constant term of the ε-expansion. 6 As it is a scalar diagram, we can make direct comparisons with FIESTA2 [33], which is much faster but applies to scalar integrals only. For AMBRE/MB the complete calculation takes about two minutes only, which is comparable to using FIESTA2.
Finally, as a more complicated case, let us take a pentabox of rank R = 3:  As it can be seen in MB_PBox.m, the construction of the MB-representation for (5.5) starts with momentum l 2 , leading to maximally nine dimensional MB-integrals. Applying Barnes lemmas reduces the integrals to at most seven-dimensional. Starting with the integration over l 1 , we get maximally eleven-dimensional integrals. Altogether, the calculation takes about five hours: The generation of the MB-representations needs a couple of minutes, and the rest of the time is needed for analytic continuations and numerical integrations.
Additional instructive examples for using AMBRE and CSectors may be found at the webpages [20,32].

Summary
The CSectors package has been prepared as an interface to the sector decomposition package sector_decomposition for the automatic numerical evaluation of tensor Feynman integrals. The AMBRE package has been extended for the automatic treatment of tensor structures of multiloop Feynman integrals.
For CSectors, the bottleneck is the numerical evaluation based on sector_decomposition and Ginac which, especially for higher rank tensors and higher loop orders, consumes a huge amount of RAM. However, for smaller problems, up to two loops and tensors of moderate rank, the program works well. For more complicated cases we recommend to use the Mellin-Barnes approach. Here there are many possibilities to optimize the way to get numerical results. This can be done at the level of construction of MB-representations (e.g. by a change of order of the integration over internal momenta), then there are different ways of analytic continuation, and finally the Mathematica package barnesroutines [34] can be used to try to reduce the dimensionality of the MB-integrals; examples for the latter can be found in [35] and at [20,32]).
For the near future it is planned to automatize the construction of MB-representations for non-planar Feynman integrals, what deserves to leave the loop-by-loop approach. Further, it is foreseen to build MB-representations for special forms of linear propagators which appear e.g. in heavy quark effective theory (HQET) or in calculations of the QCD static potential.

Appendix
The instructive examples for using AMBRE and CSectors may be found at the webpages [20,32].
invariants must be defined before Fullintegral.
This version of AMBRE uses a semi-automatic approach when building Mellin-Barnes representation. That methodology is accomplished by the following two functions.
(3) IntPart[iteration, options] -prepares a subintegral for a given internal momentum by collecting the related numerator, propagators and integration momentum: • iteration: iteration for which subintegral will be prepared. In practice IntPart function must be executed in specific order i.e firstly IntPart [1] then IntPart [2] and so.
• options: -Text: it can have two boolean values True or False. Controls if additional text appears during calculations.
(4) SubLoop[integral, options] -determines for the selected subintegral the U and F polynomials and an M-B representation.
• integral: this argument must be left as it is.
• options: -Text: as in the IntPart function.
-Xintegration: controls whether integration over Feynman parameters is performed or not.
• value: can be 0 or 1. For the first one user can modify F polynomial. For the latter this possibility is turned off.
(6) BarnesLemma[representation,number,options] -function tries to apply Barnes's first or second lemma to a given representation.
• representation: M-B representation to be checked.
• number: it has two possible values, 1 for the first Barnes lemma and 2 for the second lemma.
• options: there are the following two boolean options for this function -Text: displays or does not display an additional information.
-Shifts: searches for the pairs of two integration variables z i + z j and z i − z j which, after application of the appropriate shift one of it is cancelled. • options: -Text: displays or not information text.
-OptimizedResult: the final MB representation is written in such a way that Gamma functions are factorized. This option makes an output more condense. However sum of different Gamma functions for a given integral can cause problems if treated as it stands in analytic continuation (finding rules). Optionally it is switched off.
-BarnesLemma1: turns on or off first Barnes lemma checking.
-BarnesLemma2: does the same as above option but for second Barnes lemma Intermediate output of a given subloop is displayed in a specific format using INT function: INT[numerator,representation,propagators1,propagators2] • arguments: numerator: it is, for tensor integrals, of the form {k2[mu1],k2[mu2]}.
representation: Mellin-Barnes representation obtained during the calculation of the current sub-loop/iteration part. It is multiplied by Mellin-Barnes representations which were calculated in the previous step/iteration.
-propagators1: propagators extracted out of the F-polynomial in the current iteration.
-propagators2: keeps propagators in the same form as propagators1 does. The propagators2 are independent of the present integration variable but undergo next iteration(s).
(2) BarnesLemma[representation,number,powers of propagators,options] -this function works as in the previous versions of AMBRE. The only difference is the extra parameter: • powers of propagators: -here a list of powers of propagators which appear in an input loop integral must be given, e.g. {n1,n2,n3}.

CSectors
The appropriate loop integral is defined by: (1) a list of kinematic invariants, e.g. invariants = {p1*p1 -> s}, (2) DoSectors[numerator, propagators, internal momenta, options][min, max] The arguments of this function are exactly the same as in case of MBrepr in the AMBREnLOOP package. However, in contrast to the AMBRE package, numerator can be written only in the uncontracted form.
The package allows to modify its behaviour by adding additional options: • SetStrategy: chooses one of the strategies available in sector decomposition libraries [23].
• SourceName: a prefixing for source, binary and log files; the option just allows to choose any name for the files connected with the calculation of a given integral.
• TempFileDelete: by default it is set to TempFileDelete->True; when set to False, it does not delete C++ source and binary files as well as log file.
• LogFile: causes CSectors to create (or not) a log file, where numerical results for given integral and epsilon are stored. The default is True • ShowErrors: controls whether the errors of the numerical calculation will be displayed or not; errors are calculated using the function res.get_error() of sector decomposition [23].
• compiler: allows to choose another compiler; the default is g++.
• cppflags, libs: the paths to header and library files, required by sector decomposition libraries [23].
• min, max: indicates minimum and maximum of the Laurent series expansion in epsilon.
The default options can also be displayed by the command Options[DoSectors].