New prospects for the numerical calculation of Mellin-Barnes integrals in Minkowskian kinematics

TL;DR

This paper discusses recent advances and new prospects in the numerical calculation of Mellin-Barnes integrals in Minkowskian kinematics, highlighting methods, applications, and future potential for multi-loop Feynman diagram computations.

Contribution

It introduces a powerful automated approach for numerical MB integral evaluation using the MB-suite, and explores alternative contour deformation methods in Minkowskian kinematics.

Findings

Validated the MB method against sector decomposition in complex physical cases.

Developed a top-bottom automated approach using the MB-suite for multidimensional integrals.

Discussed the potential of MB methods for future high-energy physics applications.

Abstract

During the last several years remarkable progress has been made in numerical calculations of dimensionally regulated multi-loop Feynman diagrams using Mellin-Barnes (MB) representations. The bottlenecks were non-planar diagrams and Minkowskian kinematics. The method has been proved to work in highly non-trivial physical application (two-loop electroweak bosonic corrections to the $Z \to b \bar{{b}}$ decay), and cross-checked with the sector decomposition (SD) approach. In fact, both approaches have their pros and cons. In calculation of multidimensional integrals, depending on masses and scales involved, they are complementary. A powerful top-bottom approach to the numerical integration of multidimensional MB integrals is automatized in the MB-suite AMBRE/MB/ MBtools/MBnumerics/CUBA. Key elements are a dedicated use of the Cheng-Wu theorem for non-planar topologies and of shifts and deformations of the integration contours. An alternative bottom-up approach starting with complex 1-dimensional MB-integrals, based on the exploration of steepest descent integration contours in Minkowskian kinematics, is also discussed. Short and long term prospects of the MB-method for multi-loop applications to LHC- and LC-physics are discussed.

Figures

Tables

Table 1. Table 1: Performance tests of MBDE for the integrand F 1 ​ ( z ) = ( − s ) − z ​ Γ 3 ​ ( − z ) ​ Γ ​ ( z + 1 ) / Γ ​ ( − 2 ​ z ) subscript 𝐹 1 𝑧 superscript 𝑠 𝑧 superscript Γ 3 𝑧 Γ 𝑧 1 Γ 2 𝑧 F_{1}(z)=\left(-s\right)^{-z}\Gamma^{3}(-z)\Gamma(z+1)/\Gamma(-2z) . The relative error δ an subscript 𝛿 an \delta_{\mathrm{an}} is defined as δ an = | ( I an − I MBDE ) / I an | subscript 𝛿 an subscript 𝐼 an subscript 𝐼 MBDE subscript 𝐼 an \delta_{\mathrm{an}}=|(I_{\mathrm{an}}-I_{\texttt{MBDE}})/I_{\texttt{an}}| . I an subscript 𝐼 an I_{\mathrm{an}} is the analytical value of the integral I 1 ​ ( s ) subscript 𝐼 1 𝑠 I_{1}(s) , I MBDE subscript 𝐼 MBDE I_{\texttt{MBDE}} is the numerical value of this integral evaluated with the MBDE method. Finally, T MBDE , MB , MBnum subscript 𝑇 MBDE MB MBnum T_{\texttt{MBDE},\,\texttt{MB},\,\texttt{MBnum}} display runtimes (in seconds) needed to numerically evaluate an integral using MBDE , MB.m and MBnumerics.m (with default settings), respectively.

$s$	$I_{\text{MBDE}}$	$-{\log }_{10}{\delta }_{\mathrm{an}}$	$T_{\text{MBDE}}\left[\text{s}\right]$	$T_{\text{MB}}\left[\text{s}\right]$	$T_{\text{MBnum}}\left[\text{s}\right]$
$-1/20$	$4.96\times {10}^{-2}$	$6/9/11$	1.31/2.48/16.94	1.43	1.16
$1+i{0}^{\pm }$	-1.21	6/11	15.05/53.19	–	1.28
$5$	$4.30+14.05i$	10	13.5	–	1.57

Equations

$$\frac{1}{(A_1+\dots +A_n{)}^{\lambda }}$$

$$I(s)=\frac{1}{2\pi i}\underset{{\mathcal{C}}_0}{\int }dzF(z)=\frac{1}{2\pi i}\underset{c_0-i\infty }{\overset{c_0+i\infty }{\int }}dz{e}^{-f(z)}.$$

$${\int }_{{\mathcal{C}}_0}dz{e}^{-f}=\stackrel{\mathrm{Overall}\mathrm{factor}}{\overbrace{e^{-i\mathrm{Im}f{|}_{\mathcal{J}(z_{\ast })}}}}\stackrel{\mathrm{Damping}\mathrm{factor}}{\overbrace{\underset{\mathcal{J}(z_{\ast })}{\int }dz{e}^{-\mathrm{Re}f}}}+\stackrel{\mathrm{Remnants}}{\overbrace{2\pi i\sum _{{\mathcal{C}}_0\to \mathcal{J}(z_{\ast })}\mathrm{Res}{e}^{-f}}}.$$

$$\frac{d{z}_i(t)}{dt}=-{\left(\frac{\partial f(z)}{\partial z^i}\right)}^{\ast },z(+\infty )=z_{\ast }.$$

$$\frac{dx(t)}{dt}=\frac{\partial \mathrm{Im}f}{\partial y},\frac{dy(t)}{dt}=-\frac{\partial \mathrm{Im}f}{\partial x}.$$

Full text

\PACS

02.70.Wz, 12.15.Lk,12.38.Bx

1 Introduction

Historically a concept of Feynman diagrams was presented for the first time at a special by-invitation-only meeting at the Pocono Manor Inn in Pennsylvania in 1948 by Feynman as an alternative to procedures of perturbative calculations in QED [1, 2]. The idea was systematically treated for the first time by Dyson in his two seminal papers [3, 4] followed by Feynman himself [5, 6]111That is why initially it was being called the Feynman-Dyson approach to QED.. Integrals which stand behind the diagrams are, together with a renormalization procedure, in the core of the technical difficulties, which increase with the number of ”legs” and ”loops” involved in calculation of contemporary QCD and electroweak processes. It is clear that steady progress in particle physics needs new ideas and crafting ever-changing theoretical tools and techniques of calculations.

In the following the MB-suite will be described to some detail. It comprises several tools for dimensionally regulated Feynman integrals in the momentum space: (i) Transform them into Feynman integrals expressed by Feynman parameters (textbook knowledge); (ii) Use AMBRE [7, 8, 9, 10] — transform them into Mellin-Barns integrals, valid at initial parameters which include a finite shift $\epsilon$ of dimension, $d=4-2\epsilon$, and with original integration paths parallel to the imaginary axis; (iii) Use MB.m or MBresolve.m [11, 12] — perform an analytical continuation in $\epsilon$, approaching small $\epsilon$ and (iv) — expand the Mellin-Barnes integrals as series in small $\epsilon$; (v) Use barnesroutines [13] — perform simplifications using Barnes lemmas. (vi) At this stage the original representation of the Feynman integral in terms of several finite MB-integrals has been formulated. One may now start to calculate them, either analytically or numerically, or in a mixed approach. In sufficiently complicated situations, only numerics can be applied. (vii) Use MBnumerics.m [14] to perform parametric integrations of the MB-integrals along the paths defined in step (iii), thereby applying a variety of techniques: integration variable transformations, reparameterizations, contour deformations, contour shifts and whatsoever. For the parameter integrations, CUHRE of the package CUBA [15] is used. To some extent, we gave some descriptions of details before [9, 16, 17].

In this article we focus on purely numerical approaches to Feynman integrals developed in last few years beyond one-loop (NLO) perturbation. One is faced with several technical obstacles. There are infrared singularities. We know about two methods to treat them properly without limitations. One is the MB-method, the other one sector decomposition [18, 19], which is also numerical. We aim at direct calculations in Minkowskian kinematics, which presents serious convergence problems but is crucial for production processes at high energy accelerators like LHC and LC. No doubt that framed with powerful fast, stable, accurate and universal software, direct numerical calculations will become necessary for practical applications on mass scale, similarly as it happened at the NLO level222To solve the integrals, analytical methods can be used, though they exhibit natural limitations when sophisticated integrals with many parameters appear. Such a situation takes place in gauge theories, like in the electroweak-QCD Standard Model. However, concerning analytical approaches to Feynman integrals, we should especially appreciate recent progress in differential equation method [20, 21, 22, 23], which got a push in 2013 [24] followed by latest corresponding software and ideas [25, 26, 27, 28, 29]. Here further progress in developing integration-by-parts (IBP) concepts is also very important [30]. , see FeynArts/FormCalc [31, 32], CutTools [33], Blackhat [34], Helac-1loop [35], NGluon [36], Samurai [37], Madloop [38], GoSam [39], PJFry [40], OpenLoops [41] and [42, 43, 44, 45].

In this article we discuss the state of the art of purely numerical approaches to multiloop integrals, focusing on the Mellin-Barnes method. We show the first completed and non-trivial application in cutting-edge physical calculations using the MB-suite followed by further perspectives.

2 Numerical concepts beyond NLO level

Fully numerical techniques for the evaluation of two- and higher-loop integrals need the extraction of ultraviolet, infrared and collinear singularities. On top of that, they must be numerically stable and efficient. A qualitative comparison of different numerical integration techniques for Feynman parameter integration of massive multi-loop integrals can be found in [46, 47, 48, 49]. The main methods are dispersion relations, Bernstein-Tkachov method, differential equations, subtraction terms, the SD and MB methods. Here we will discuss the last two, specifically focusing on MB. There are presently only few public programs for the numerical integration of integrals beyond NLO level. NICODEMOS [50] is based on contour deformations. There are also complete programs dedicated specifically to the precise calculation of two-loop self-energy diagrams [51, 52]. However, the most advanced and universal programs are based on the SD or MB approaches: Sector decomposition, developed into independent packages (present versions) Fiesta 4 [53] and SecDec 3 [19] followed by pySecDec [54]; MB.m [11] and MBresolve.m [12] packages extract $\epsilon$-singularities in dimensional regularization of MB multiloop integrals, offering also possibilities of numerical integrations in Euclidean kinematics, which is relatively simple as no physical branch cuts are present there. It was used intensively in the past to cross check numerically analytical results for multiloop integrals. In next section we will discuss new ideas for making possible MBnumerical integrations directly in the physical region.

3 Numerical integrations of MB integrals in Minkowskian region

The Mellin-Barnes transformation of Feynman multidimensional integrals to multivariable complex plane integrations [55, 56] has been used in many particle physics calculations. In the first applications [57, 58] this kind of transformation has been applied directly to propagators in the loop integrals, changing $"\rm{momenta}^{2}-\rm{mass}^{2}"$ terms into ratios of momenta and masses in complex plane. Nowadays, a more efficient and systematic treatment of multiloop integrals goes by expressing Feynman integrals by the Symanzik polynomials $F$ and $U$ [59, 60, 61], for which the general MB formula is applied

[TABLE]

As we can see, $n$ additive terms lead to $n-1$ complex integrals. The $A_{i}$ terms correspond to kinematical parameters of the integral. A typical simple example is the 1-dimensional singular part of the 1-loop massive QED vertex [11, 9, 62] $\sim\int dz(-s)^{-z}\Gamma^{3}(-z)\Gamma(1+z)\Gamma^{-1}(-2z).$ Choosing properly the contour of integration can make the annoying oscillatory behavior of the term $(-s)^{-z}$ small and controllable (for $s>0$, so Minkowskian kinematic). Furthermore, Gamma functions $\Gamma$ exhibit singularities either, and make the task of integral evaluations highly non-trivial.

The construction of MB integrals through Symanzik polynomials is automatized in the AMBRE project [7, 8, 9, 10]. Using it with MB.m or MBresolve.m, IR and UV divergencies can be extracted and regulated multidimensional MBintegrals are obtained [63]. On the webpage [13] more auxiliary packages with examples related to MB calculations can be found.

The first serious trial directed to the numerical integration of MB integrals in Minkowskian space-time was undertaken in [64]. The method developed there is based on simultaneous rotations of integration paths for all variables by the same angle in the complex plane and has been applied successfully to the calculation of two-loop diagrams with triangle fermion subloops for the $Z\to b\bar{b}$ formfactor [65]. Another interesting numerical application of MB integrals for phase space integrations can be found in [66] and [67, 68]. There some parametric integrals are considered and transformations of MB integrals into Dirac delta constraints have been explored.

Now we will present recent developments. First we describe a Top-Bottom approach in which the MBnumerics.m package deals with multidimensional MB integrals; it was described partly in [9] and applied in [17]. Another Bottom-Top approach is at the exploratory stage; optimal complex contours of MB integrations are worked out systematically for one-dimensional MB integrals [62].

3.1 Top-Bottom approach - shifts, deformations and MBnumerics.m

As we can see from (1), Gamma functions are there [56]. In Fig. 1 the real part of the $\Gamma[z]$ function is sketched. It is regular in positive $\Re[z]$ and has singularities for integer negative $\Re[z]$.

Note also that at the negative axis between the pole positions, the integrand becomes smaller in its absolute value for the function evaluated at an argument further away from the origin. In addition, for a pole crossed by an argument shift, one has to add the corresponding residue which by itself is also an integral, but will have a dimension less than the original one. Repeating the procedure for several integration variables, the original MB integral gets replaced by several lower-dimensional integrals, and the original one with a shifted integration path. In the end, the (module of the) resulting contribution of the original integral after shifts can be made smaller than the desired accuracy of the calculation. In effect, the procedure constructs a sum over a finite number of residues with a controlled remainder. This procedure of shifts is implemented in MBnumerics.m [14]. Some other important features of the procedure like contour deformations and mappings of integrated variables into finite intervals have been discussed in [9, 16]. Fig.2 sketches roughly the idea.

In the project [17] we derived Mellin-Barnes representations for all integrals, which had up to eight dimensions. For a cross-check, each integral was computed with MB and SD techniques. There are only few classes of diagrams for which eight digits could not be achieved with both methods, an example is given in Fig.3; for further discussion, see [9, 16]. These diagrams have high order divergiences and an application of the sector decomposition approach leads to numerical problems related to an accuracy and time consumption. In contrast to this, the corresponding MB integrals can be computed with reasonable computer time resources.

Typically, for integrals which involve many masses, SD fits better while the MB method works out perfectly for more ”massless” diagrams. Thus, the MB suite and the sector decomposition techniques are to a large extent complementary [8] and both numerical methods can be successfully explored together in cutting-edge physical problems.

In [17] for MB and planar diagrams the newest version AMBRE v2.1 [10] is used, for non-planar diagrams it is AMBRE v3.1 [10]. Planarity of diagrams is controled automatically with the PlanarityTest.m package [71, 72]. Numerical results have been obtained using MBnumerics.m [14]. As it is demonstrated in Fig. 2, the shifts accumulate at each new iteration many residues, until the desired accuracy is reached. It is worth noting that the integration error of MBnumerics.m is mostly dominated by the collection of residues which have fast convergence. For higher-dimensional integrals, MBnumerics.m collects more residues. The resulting error from all residues is determined by Pythagorean addition. In Fig. 2 the double arrows with zero marks denote pairs of residues which are identified to finally cancel exactly. To identify such pairs to a high accuracy, MBnumerics.m performs the integration of the corresponding candidates at a different kinematical point where a high numerical accuracy is reached. If then the integrals agree up to a sign, MBnumerics.m sums them up to zero. This is only one example of many numerical problems which have been solved in the MBnumerics.m algorithm, in order to get highly accurate numerical results in the Minkowskian region. The package is yet under development, and our present estimation is that in the near future even 12-dimensional MB integrals can be touched – e.g. pentaboxes.

3.2 Bottom-Top MBDE approach - optimal steepest descent contours

In a nutshell, this is a stationary phase method leading to optimal steepest descent integration contours. They can be found using Lefschetz thimbles (exact contours) or their Padé approximation [62].

Lefschetz thimbles (LT) are a fascinating subject, crossing many issues like behaviour of LT in presence of poles, singularities and branch cuts, behaviour in the complex infinity, Stokes phenomenon, relation to relative homology of a punctured Riemann sphere, etc. It can be applied e.g. to the analytical continuation of 3d Chern-Simons theory, QCD with chemical potential, resurgence theory, counting master integrals or repulsive Hubbard model. Still, applying this method to the numerical evaluation of MB integrals is at the exploratory stage and an effective and general determination of multivariate MB contours must be worked out yet in more detail.

In this section we present the main idea as an alternative approach to the numerical computation of MB integrals, starting from the bottom, the lowest one-dimensional MB integrals, in both Euclidean ($s<0$) and Minkowski $(s>0)$ regions. These cases have been explored in fine details in [62].

Let us write a general MB integrand $F(z)$, transformed into exponential form:333For brevity, we suppress the dependence on $s$ and shall use $F(z)$ instead of $F(s,z)$.

[TABLE]

$\mathcal{C}_{0}$ is a contour defined by $\mathrm{Re}\,z=c_{0}$ while $f(z)=-\ln F(z)$.

One of possible ways to get rid of numerical problems with the MB integrand $F(z)$ which is of highly-oscillatory behaviour [9] is to integrate (2) over a new contour $\mathcal{C}=\mathcal{J}_{1}+\mathcal{J}_{2}+\mathcal{A}$. A typical example is sketched in Fig. 4 where $\mathcal{C}$ is a sum of three contours $\mathcal{J}_{1}$, $\mathcal{J}_{2}$ and $\mathcal{A}$ along which the behaviour of $f$ is under control.

Taking $f=\mathrm{Re}\,f+i\,\mathrm{Im}\,f$, we deform the original integration contour $\mathcal{C}_{0}$ to a Lefschetz thimble $\mathcal{J}(z_{*})$,

[TABLE]

The analytical formula describing $\mathcal{J}_{k}$ can be found only in the simplest cases by explicit solving the equation $\mathrm{Im}\,f=\mathrm{const.}$ Instead, we use the fact that the function $\mathrm{Re}\,f$ defines a Morse flow [73, 74]. Such a flow is realized by a parameterization $t\mapsto z(t)$ of $\mathcal{J}_{k}(z_{*})$ in a form of Lefschetz thimbles [75, 76, 77, 78, 79, 80]. The Lefschetz thimble $\mathcal{J}(z_{*})$ is defined as a union of curves $t\rightarrow z(t)=(z_{1}(t),\ldots,z_{i}(t),\ldots,z_{n}(t))\in\mathbb{C}^{n}$ which satisfy the following differential equation [76, 79]:

[TABLE]

Here $z_{*}$ is a saddle point of a meromorphic function $f$. The crucial observation is that for $\mathcal{J}(z_{*})$ we can take $\mathrm{Im}\,f$ = const, leading to the overall factor in (3). Note that $\mathrm{Im}\,f$ generates a Hamiltonian flow on $\mathbb{R}^{2n}$, e.g. for $n=1$:

[TABLE]

The $\mathrm{Re}\,f$ is monotonically decreasing when $t\rightarrow+\infty$ and goes to $+\infty$ when $t\to-\infty$, leading to the damping factor in (3).

Remnants in (3) are, according to Cauchy’s theorem, residues over poles when the integration contour is deformed from $\mathcal{C}_{0}$ to $\mathcal{C}$ and encircles extra poles of $F=e^{-f}$, as shown in Fig. 4. In this way, (4) gives a possibility to solve for $z(t)$ such that the integral (3) is under control. As $\mathcal{J}_{1,2}$ we choose such stationary phase contours which start at saddle points $z_{*}^{(1,2)}$ and go towards infinity without hitting other poles. Both contours are chosen such that $\mathrm{Im}\,f$ is constant along them and function $\mathrm{Re}\,f$ is strictly increasing when one moves away from $z_{*}^{(1,2)}$. Varieties defined in such a way are called steepest descent contours [81, 82]. Usage of $\mathcal{J}_{1,2}$ allows to control the behavior of $f(z)$ when $z\to\infty$. Because $\mathrm{Re}\,f$ is stricly increasing, the integrand $e^{-f}$ decreases rapidly at the ends of $\mathcal{J}_{1,2}$. That transform the integral (2) into a form which is more suitable for numerical treatment.

With respect to various methods known in the literature [48, 11, 9, 64] which shift/rotate contours or use approximate forms thereof, the MB method which relies on the differential equation (4), in short the MBDE method, relies on deriving the numerical parameterization $z(t)$ of $\mathcal{J}_{k}$ as a solution of the differential equation (4) and then, again numerically, integrating the function $e^{-\mathrm{Re}\,f}$ along the contour $\mathcal{C}$ composed of Lefschetz thimbles $\mathcal{J}_{k}$ (and the compact contour $\mathcal{A}$ if necessary). The purely numerical approach MBDE is complementary to the Padé approximation presented in [62].

Let us shortly discuss numerical features of the MBDE metod and display results of some performance tests. We stress that the tests are preliminary, implemented directly in Mathematica, in graphical mode, on an i7 $2.9$ GHz CPU. Both kinematical regions $s<0$ and $s>0$ are treated in the same way in MBDE, although $s>0$ seems to be more CPU time consuming. For a final accuracy of the order of $10^{-6}$ the method is as fast as MB [11] and MBnumerics [83], while for an accuracy of $10^{-11}$ and higher, MBDE turns out to be more than 10 times slower than other two packages; see Tab. 1. To get a precision of the order of $10^{-16}$ some kind of optimization of the method is needed. Presumably, it can be made much faster by implementing in e.g. Fortran or C/C++ or by applying a dedicated method of solving differential equations. Parallelization or dividing integration regions into smaller parts to achieve larger precision are also possible options.

4 Summary and Outlook

In the last few years there is substantial progress in the direct calculation of multiloop integrals (Feynman diagrams) in the physical, Minkowski regime using both SD and MB methods. Both methods are complementary in several respects. In the MB case, the most advanced is the top-bottom approach implemented in the MBnumerics.m package where multidimensional MB integrals can be solved in physical kinematics with high accuracy for MB integrals of dimension eight and below. Potential applications of the discussed numerical methods are complete 2-loop electroweak pseudoobservables needed for future linear colliders – multi-massive 2-loop vertices – and also non-resonant two-loop box diagrams – complete cross sections, including LHC problems [64]).

Using numerical methods, we are approaching automation in calculation of Feynman integrals beyond the NLO level directly in physical kinematics. Perspectives are robust, concerning both high and low energy physics.

Acknowledgements

I.D. is supported by a research grant of Deutscher Akademischer Austauschdienst (DAAD) and by Deutsches Elektronensychrotron DESY; J.G. and T.J. are supported by the Polish National Science Centre (NCN) under the Grant Agreement No. DEC-2013/11/B/ST2/04023; T.R. is supported in part by an Alexander von Humboldt Polish Honorary Research Fellowship. The work of J.U. is supported by Graduiertenkolleg 1504 “Masse, Spektrum, Symmetrie“ of Deutsche Forschungsgemeinschaft (DFG).