On the Resolution of Singularities of Multiple Mellin-Barnes Integrals

One of the two existing strategies of resolving singularities of multifold Mellin-Barnes integrals in the dimensional regularization parameter, or a parameter of the analytic regularization, is formulated in a modified form. The corresponding algorithm is implemented as a Mathematica code MBresolve.m


Introduction
The method of Mellin-Barnes (MB) representation is one of the most powerful methods to evaluate multiloop Feynman integrals. It is based on a very simple formula, which is called MB representation and is applied to replace a sum of two terms raised to some power by their products in some powers. The contour of integration should be chosen so that the poles with a Γ(. . . + z) dependence are to the left of the contour and the poles with a Γ(. . . − z) dependence are to the right of it.
The simplest way to apply (1) is to represent massive propagators as continuous superpositions of massless ones, at the cost of introducing additional integrations over MB variables. More often, MB integrals are introduced at the level of Feynman or alpha parameters. Anyway, (1) is used in an appropriate way, with the goal to obtain integrals over loop momenta or Feynman/alpha parameters which can be taken in terms of gamma functions. As a result one obtains a multiple MB integral with an integrand involving gamma functions and powers of kinematic invariants.
It is important to derive a MB representation for general values of the powers of the propagators. First, this provides unambiguous prescriptions for a choice of contours of integration. The resulting rule is similar to the rule for (1): when we understand a multiple MB integral iteratively and consider an integration over a variable z i it is implied that the poles with a Γ(. . . + z i ) dependence are to the left of the contour and the poles with a Γ(. . . − z i ) dependence are to the right of it. Second, such a general MB representation provides crucial checks and can be applied at any chosen values of the powers of the propagators.
When evaluating a Feynman integral for specific powers of the propagators (indices), one starts from the general MB representation and obtains an integral of the form where ε = (4−d)/2 is the dimensional regularization parameter, a i , . . . , c ′ ij are rational numbers, x k are ratios of kinematic invariants and/or masses, and their exponents, d k , are linear combinations of ε and z-variables. Typically, c ij = ±1.
Although results for dimensionally regularized Feynman integrals are practically needed in a Laurent expansion in ε, one could try to evaluate a given MB integral (2) for general ε. Then the ε-expansion is performed in the result. However, for sufficiently complicated multiloop Feynman integrals, the evaluation at general ε turns out to be impossible, and one proceeds with an ε-expansion.
One of the advantages of this method is that the singularity structure in ε can be resolved in a simple way. This procedure basically consists of taking residues and shifting contours, with the goal to obtain a sum of integrals where one can expand integrands in Laurent series in ε. To do this one can apply two strategies formulated in [1] and [2] which we are going to call Strategy A and Strategy B, respectively.
Strategy A is described and illustrated in numerous examples in Chapter 4 of [3]. It was applied, e.g., in [4]. According to strategy A, one performs an analysis of the integrand to reveal how poles in ε arise. The guiding principle is that the product Γ(a + z)Γ(b − z), where a and b can depend on the rest of the integration variables, generates, due to the integration over z, the singularity of the type Γ(a + b). Indeed, if we shift an initial contour of integration over z across the point z = −a we obtain an integral over a new contour which is not singular at a + b = 0, while the corresponding residue involves an explicit factor Γ(a + b). This observation shows that any contour of one of the following integrations over the rest of the MB variables should be chosen according to this dependence, Γ(a + b). Hence one thinks of integrations in various orders and then identifies some 'key' gamma functions which are crucial for the generation of poles in ε. Then one takes residues and shifts contours, starting from the first poles of these key gamma functions. The same analysis and procedure is applied to the contributions of the residues.
Within Strategy B, one chooses an initial value of ε and values of the real parts of the integration variables, z i , . . . in such a way that the real parts of all the arguments of the gamma functions in the numerator are positive and one can integrate over straight lines. Then one tends ε to zero and whenever the real part of the argument of some gamma function vanishes one crosses this pole and adds a corresponding residue which has one integration less and is treated as the initial integral within the same procedure.
Strategy B is algorithmic in its character and, indeed, two algorithmic descriptions were formulated in [5,6]. A public code called MB.m was presented in [6]. Strategy B was successfully applied, e.g., in [7,8] and in many other papers.
The purpose of this letter is to formulate Strategy A in a slightly modified form and to present the corresponding algorithm implemented in Mathematica.

The modified Strategy A
To present a modified Strategy A let us explicitly formulate what was implied in the initial Strategy A. When we take care of one of the key gamma functions we shift a contour and take a residue. Let Γ(A i ) with A i = a i + b i ε + j c ij z j be one of the key gamma functions in (2). Without loss of generality, we can consider ε real. Then changing the nature of the first pole of this gamma function means changing the rule for an admissible contour, i.e. that, instead of the condition ReA i > 0 when crossing the real axis in the process of the integration, we have the condition −1 < Re A i < 0. Let us denote this transition by replacing Γ(A i ) by Γ (1) (A i ). The initial rule for the contour can be changed again and then we have the condition −n < Re A i < −n + 1 for n = 2, 3, . . . with the notation Γ (n) (A i ).
Within Strategy B one has straight contours in the beginning. Rather, in the modified Strategy A, we will be oriented at straight contours in the end. Apparently, it is desirable to achieve a minimal number of terms after the resolution of the singularities in ε. To do this, let us try to search for contours which are going to have in the end of this procedure and for which the gamma functions in the numerator are changed, in the above sense, in a minimal way.
To formalize this requirement, let us introduce the function is the integer part of a number and x + = x for x > 0 and 0, otherwise. In other words, if −n < x < −n + 1 then σ(x) = n for n > 0 and σ(x) = 0 for n ≤ 0.
So, let us set ε = 0 and search for contours, i.e. Rez i , for which the sum After such a choice is done we identify gamma functions which should be changed, in the above sense, in order to arrive at a final integral where a Laurent expansion in ε is possible. In fact, this step replaces the first step in the primary Strategy A where one identified such key gamma functions after the analysis characterized above.
Then the second step in Strategy A is the same as in the old version: we take care of the distinguished gamma functions, i.e. take a residue and replace Γ by Γ (1) (A i ) (and, possibly, Γ (1) (A i ) by Γ (2) (A i ) etc.) We proceed iteratively, as in the previous strategy: every residue is considered from scratch, i.e. treated in the same way as the initial MB integral.
Let us emphasize that although this strategy aims to minimize the number of resulting terms, we cannot exclude that there is another way of resolving the singularities in ε that is the best one in this sense. For example, it can happen, in rather complicated examples, that different orders of changing the key gamma functions lead to different numbers of resulting terms. Still we believe that such a difference is negligible and that the Strategy A provides a resolution of the singularities in ε at least very close to the theoretically best one.
The difference of the new and the old Strategies A is minor.
In the examples of [4] one can see that resulting contours were straight indeed. This difference can still be seen in the following simple example of the integral which, of course, can be evaluated at general ε by the first Barnes lemma. Within the 'old' Strategy A, one can observe that there is no gluing of poles so that one can expand the integrand in ε. However, a resulting contour cannot be chosen as a straight line. Rather, within the 'new' Strategy A, we have to choose a gamma function to be modified and, as a result, we obtain a residue and an integral over a straight line which both can be expanded in ε.

Implementation in Mathematica
The code should be loaded together with the package MB [6] by Czakon. It uses the MBresidue routine and some other private functions from that package; after obtaining the list of integrals one can continue the contour optimization and the numerical evaluation with the standard routines of MB.
Formally the code can be described the following way: the initial function results in contour prescriptions (the arguments of Γ-functions have to be positive in the end On exception move to the next For cycle 7. represent x as cα + r, where c = 0, α is one of the integration variables and r does not depend on α 8. Return t1  The function should be a function of integration variables such that the poles are only due to the Gamma and PolyGamma factors (with linear arguments). The contour prescriptions cont pr are a list of pairs x, n where x is a linear function of integration variables and n is an non-negative integer. If n is equal to zero, such a term means that x has to be positive; if n is positive, then x has to be greater that −n and smaller than −n + 1. If the cont pr parameter is missing, it is created automatically by considering all arguments of Gamma functions and their derivatives and assuming them all to be positive. An uncaught exception at the top level appears if there is a degenerate case and one is required to introduce extra regularization parameters.
Hence, instead of using The corresponding Mathematica code MBresolve.m is public and can be found at http://projects.hepforge.org/mbtools/ together with other tools for evaluating MB integrals.
Here is an example of a tenfold MB representation derived loop by loop for the four-loop ladder massless on-shell diagram with p 2 i = 0, i = 1, 2, 3, 4, where s = (p 1 + p 2 ) 2 and t = (p 1 + p 3 ) 2 .
In [

Discussion and perspectives
Let us remind that, sometimes, a given MB representation can be ill-defined even for a well-defined Feynman integral, in the sense that a gluing of poles is present, i.e. one can distinguish a subproduct of the gamma functions such that the sum of there arguments, at general ε, is equal to a non-positive integer number. To cure such a MB representation, one can introduce an auxiliary analytic regularization into some index, i.e. a i → a i + y and then analytically continue the given MB integral, first, in y to the point y = 0 and then ε to the point ε = 0. Let us stress that in such situation one starts, within the new Strategy A, with setting y = 0, ε = 0 and then searches for appropriate contours for which the gamma functions in the numerators are changed in the minimal way, similarly to the case without such analytic regularization.
Let us emphasize that both Strategy A and Strategy B are based on the fact that singularities in ε can be generated, because of gluing of poles of different nature, by the integration at compact regions. At least for planar diagrams, for which one can apply the code called AMBRE [9] based on the loop-by-loop strategy, this looks to be the only source of the singularities. However, for MB representations derived within the loop-by-loop strategy for nonplanar diagrams, there can be another source of poles. 3 This feature can be exemplified by the massless two-loop nonplanar diagram with two external legs on-shell. 4 In the corresponding MB representation derived within the loop-by-loop strategy one meets, in particular, the following onefold MB integral There is no gluing of poles so that a pole in ε cannot be generated by the integration over finite regions. Still a pole is generated and this can be seen by an explicit evaluation of this integral by closing the integration contour to the right and summing up the resulting series. This can be seen also by analyzing the asymptotic behaviour of the integrand at infinity. Setting z = x + iy and using formulae of the asymptotic behaviour of gamma functions at large arguments in the complex plane, we can observe that the leading asymptotic behaviour when y → +∞ is 1/y 1−2ε which explains the appearance of the pole. Let us still mention that, for this concrete diagram, there is a better way to obtain a 'good' MB representation -see, e.g., Chapter 4 of [3]. For simple MB integrals, say, up to six-fold, the codes MB.m and MBresolve.m work more or less in the same way. For higher-fold MB integrals, e.g. ten-fold ones, MBresolve.m produces at first fewer integrals after the resolution of singularities in ε. However, the command MBmerge of MB.m combines many terms, and the number of resulting integrals becomes smaller that produced by MBresolve.m. We should bear in mind, though, that in the current version of MB.m, this command combines all integrals which have the same integration contours, while it is more natural to consider separately integrals with different patterns of gamma functions, both from analytical and numerical points of view. Presumably, the MBmerge command can be improved in this respect.