Functional reduction of Feynman integrals

A method for reducing Feynman integrals, depending on several kinematic variables and masses, to a combination of integrals with fewer variables is proposed. The method is based on iterative application of functional equations proposed by the author. The reduction of the one-loop scalar triangle and box integrals with massless internal propagators to simpler integrals is described in detail. The triangle integral depending on three variables is represented as a sum over three integrals depending on two variables. By solving the dimensional recurrence relations for these integrals, an analytic expression in terms of the $_2F_1$ Gauss hypergeometric function and the logarithmic function was derived. By using the functional equations, the one-loop box integral with massless internal propagators, which depends on six kinematic variables, was expressed as a sum of 12 terms. These terms are proportional to the same integral depending only on three variables different for each term. For this integral with three variables, an analytic result in terms of the $F_1$ Appell and $_2F_1$ Gauss hypergeometric functions was derived by solving the recurrence relation with respect to the spacetime dimension $d$. The reduction equations for the box integral with some kinematic variables equal to zero are considered.


Introduction
Theoretical predictions for experiments at the LHC [1,2] as well as at future colliders such as the FCC [3] demand knowledge of precise radiative corrections. Precise experimental measurements are to be interpreted with sufficient precision of theoretical predictions. The problem of calculating such radiative corrections is associated, in particular, with the need to compute Feynman integrals depending on several kinematic variables and/or masses. Over the past few decades enormous progress has been made in solving the problem of evaluating Feynman integrals.
However, for further progress in the analytic evaluation of Feynman integrals, especially integrals depending on several kinematic variables and masses, new mathematical methods need to be elaborated. In this respect, more and more attention is paid to methods based on solving different kinds of recurrence relations. The first result obtained in this approach was analytic evaluation of the two-loop propagator integral with massless propagators [4]. In this paper, an analytic result for the integral was derived by solving the recurrence relation with respect to the power of a propagator. A systematic approach for evaluating Feynman integrals by solving the recurrence relation with respect to the power of a propagator was described in Ref. [5].
In Ref. [6], a method for evaluating Feynman integrals based on the recurrence relations with respect to the space-time dimension d was suggested. It turns out that hypergeometric functions appearing in the solution of the recurrence relations with respect to d [7], [8] have fewer arguments than in the results for these integrals obtained by other methods. For example, the results obtained in Refs. [9], [10] using the Mellin-Barnes integration technique, and in Ref. [11] using the negative dimension method are expressed in terms of hypergeometric functions with more arguments than those obtained for these integrals by solving dimensional recurrence relations.
In Ref. [12], new relationships between the Feynman integrals with different kinematic variables were discovered. A method of deriving functional equations from algebraic relations for products of propagators was recently proposed in Ref. [13]. At the one-loop level some functional relationships were also considered in Refs. [14], [15].
It was shown that these relationships, or in other words functional equations, can be used to express Feynman integrals in terms of integrals with fewer variables. In Ref. [16], the functional equations were used to obtain relations between integrals appearing in radiative corrections for different physical processes. An important step in evaluating radiative corrections for physical processes is the Laurent expansion in the ε = (4 − d)/2 of the analytic results for Feynman integrals. Quite essential progress in this field was made in many papers. See, for example, [17][18][19][20][21][22]. Up to now the Laurent expansion of Feynman integrals in ε is not a completely solved problem. Even at the one-loop level only the first several terms in the ε expansion of the four-and higher point functions are known. The existing results [23], [24] are not so easy to generalize for integrals depending on several masses and/or several external off-shell momenta.
It is evident that for the ε expansion the simplicity of the analytical results for dimensionally regularized integrals is rather important. For this reason the method based on the solution of recurrence relations and the method of functional reduction suggested in this article can be very useful.
In the present paper, we propose a framework for systematically reducing Feynman integrals depending on several kinematic variables and masses to a combination of integrals with fewer variables. This framework is based on solving the functional equations for Feynman integrals proposed in Refs. [12,13,25]. The main steps of our approach will be illustrated on the one-loop integrals with massless propagators.
In a sense, the application of functional equations for evaluating integrals is analogous to the use of recurrence relations with respect to some discretely changing parameters, like space-time dimension d or power of a propagator. Applying such recurrence relations, one can reduce an integral to a set of basis integrals which are in fact boundary values of the integrals of interest. Using functional equations one can reduce an integral to a combination of integrals with fewer variables, i.e. integrals defined on some hypersurfaces. In other words, these integrals can be interpreted as a kind of boundary integrals. This paper is organized as follows.
In section 2, we briefly discuss the method of discovering functional equations for Feynman integrals and describe the methods of obtaining their solutions. As an illustrative example, we consider the solution of the functional equation for the one-loop propagator integral with arbitrary masses.
In section 3, the one-loop integral associated with the triangle Feynman diagram with massless internal propagators is considered. We present the functional equation for this integral and describe its solution. The analytic result for the integral appearing in the solution of thefunctional equation is obtained as a solution of the dimensional recurrence relation. A particular case of the functional equation for the triangle integral is considered.
In section 4, we present the functional equations for the one-loop scalar integral associated with the Feynman diagram with four external legs. A two step procedure, based on functional equations, for reducing the integral depending on six variables to a combination of integrals depending on three variables is described. For these integrals, depending on three variables, an analytic result as a solution of the dimensional recurrence relation is presented. Functional reduction of the box integral for several particular cases of kinematic variables is considered. The first few terms in the Laurent expansion around d = 4 and d = 6 for these integrals are given.
In section 5, we report our conclusions and discuss future applications of functional equations for evaluating Feynman integrals corresponding to diagrams with massive internal lines and with more external legs and loops.
In Appendix A, we present definitions and explicit formulae for the Gram de-terminants and polynomials occurring in the paper. In appendix B, a derivation of the analytic result for the one-loop integral with massless internal propagators with particular emphasis on its dependence on the small imaginary part needed for the correct analytic continuation of the integral is presented. In Appendix C, the series and integral representations for the hypergeometric functions used in the paper are given.

Functional equations and their solution
At the present time, there are three methods for deriving functional equations for Feynman integrals. The method proposed in Ref. [12] is based on exploiting recurrence relations obtained by the method of generalized recurrence relations [6]. By choosing some kinematic variables, one can eliminate most complicated integrals from the recurrence relation so that the sum of remaining terms represents the functional equation. The second method is based on algebraic relations for a sum of products of propagators [13]. Integrating such sums with respect to a common to all propagators momentum gives a functional equation. The third method is based on the use of algebraic relations for modified propagators [13]. Integrating an algebraic relation depending on modified propagators with respect to a common to all propagators momentum, transforming the resulting integrals to integrals over Schwinger parameters and then mapping these integrals to the required Feynman integrals by choosing auxiliary parameters from deformed propagators lead to a functional equation.
The following questions arise naturally: how to solve the functional equations and how to use them for simplifying evaluation of Feynman integrals? We shall try to answer these questions in the next sections of this paper.

Definitions and methods of solution
A functional equation can be considered as an equation involving independent variables, known functions, unknown functions and some constants [26]. In a functional equation the unknown is a function. Rather often, the functional equation connects a function with its value for some other arguments. There is a vast literature on functional equations [26][27][28][29][30][31][32]. Solution of a functional equation is a rather difficult problem. However, there is a number of the most frequently used methods for its solution. A systematic description of such methods is given in Ref. [26]. Many methods described in this book and also in Refs. [27], [32]

Mixed methods
All these methods to some extent can be used for solving functional equations for Feynman integrals. In the present paper, the methods 1, 3 and 5 will be exploited.

Solution of the functional equation for the propagator integral
As an illustration of the first method from the above list, we shall consider the solution of the functional equation for the one-loop scalar propagator integral: where iη is the small imaginary part which fixes the analytic continuation of the integral. In Refs. [12], [13], the following relationship for this integral was derived:  2 (m 2 i , m 2 j ; s ij ), which is considered as a function of three continuous variables s ij , m 2 i , m 2 j . To solve equation (2.2), we will exploit a method, which was used for the solution of Sincov's equation [33], [34] : Setting in this equation z = 0 and assuming that the function f (x, z) is not singular at this point, we obtain the general solution Thus, using the fact that the left-hand side of equation (2.6) does not depend on z, we express the function f (x, y) as a combination of its "boundary values".
It is easy to see that the functional equation (2.2) is rather similar to Sincov's equation (2.6). Since at m 2 0 = 0 the invariants s i0 , s j0 and the integral I Therefore, relation (2.9) represents the integral depending on three variables in terms of integrals depending on two variables.
Expression (2.9) is a solution of equation (2.2) for arbitrary value of the mass m 2 0 . Indeed, substituting (2.9) in both sides of equation (2.2), simplifying arguments (2.11) of integrals, after simple algebraic transformations, we find that on the righthand side the integrals I Notice that to reduce Feynman integrals to simpler ones, the question whether expression (2.9) is a general solution of the functional equation (2.2) or not is not relevant. For our purposes it is enough to have a particular solution reducing complicated integral to a combination of integrals with fewer variables. Other sets of particular solutions will lead to another representation of the complicated integral in terms of simpler ones. These sets of integrals may be related, for example, by analytic continuation or some transformation analogous to the known transformations for hypergeometric functions.
It should be noted that x j , s j0 in equations (2.3), (2.4) and x j , s j0 in equations (2.10),(2.11) do not depend on iη that can lead to ambiguity in choosing the sign of the square root. Nevertheless, the functional equation will be valid for any choice of sign. The signs in x j , x j are to be properly correlated with the sings in s j0 , s j0 . The possibility to choose different signs of the square root means that there are two different representations of the integral in terms of simpler integrals. The integrals on the right-hand side of these two different representations depend on different sets of arguments. Excluding the initial integral from these equations will give a functional equation for integrals with fewer arguments. As it was already shown in Ref. [12], this kind of functional relations may be used for the analytic continuation of integrals with fewer variables.
We conclude this section with various remarks. First, it is interesting to note that the position of the threshold s ij = (m i + m j ) 2 for the integral on the left-hand side of equation (2.9) corresponds to the positions of thresholds s i0 = m 2 i and s j0 = m 2 j for integrals on the right-hand side.
Second, we notice that the functional equations can be used not only for reducing complicated integrals to their "boundary integrals" but also for analytic continuation of these "boundary integrals".
Third, to find an analytic expression for the simple integrals, which cannot be simplified anymore by using functional equations, one should use other computational methods. In the next sections, we will use dimensional recurrence relations for the triangle and box integrals at the final stage of calculation.

Functional reduction of the integral I (d) 3
In the present paper, the functional equation and its solution for the one-loop triangle integral with all internal masses equal to zero will be considered. To solve the functional equation for this integral, we will use the functional equation for the integral with massive lines: where 3.1 Derivation of functional equation for the integral I (d) 3 and its solution To derive the functional equation for the integral (3.1), one can exploit the algebraic relation for the products of three propagators [13]: As it was shown in Ref. [13], relation (3.3) is valid if 4) and the parameters m 2 0 , x j obey the following system of equations: Solving this system of equations for x 1 ,x 2 , we have where Λ 3 is the root of the equation Integration of relationship (3.3) with respect to the momentum k 1 gives a functional equation for the one-loop integral I (d) 3 with arbitrary masses: where s 12 , s 13 , s 23 are independent scalar invariants and s j0 are given by Setting m 1 = m 2 = m 3 = m 0 = 0 in equations (3.10), (3.11) and replacing s ij → q ij in order to avoid confusion in the notation, give the functional equation for the massless case 3 (0, 0, 0; q 20 , q 10 , q 12 ), (3.12) where q 10 = q 12 λ 3 + q 13 z 3 , the parameter z 3 is arbitrary and λ 3 is a solution of the quadratic equation: 14) The right-hand side of equation ( by choosing z 3 . To find a solution of equation (3.12), one can use an approach described in [26]. Namely, we will find a solution for the integral with massless propagators from a more general functional equation which can be obtained from formula (3.10). As was noted in Ref. [32], it can happen that solving a more general equation can be easier than solving a particular case of this equation. Setting m 2 1 = m 2 2 = m 2 3 = 0 in equation (3.10) but keeping m 0 different from zero, we find: where and Λ 3 is the solution of the equation: We note that in the obtained equation there are more arbitrary parameters than in equation (3.12) and also new kinds of integrals, namely, integrals with one massive internal line. It is expected that the solution sought may be found by choosing arbitrary parameters x 3 , m 0 . Compared with Sincov's equation, it is not so easy to find values of the arbitrary parameters leading to the reduction in the number of independent variables for all the integrals simultaneously. The reduction in the number of variables can take place if by choosing parameters x 3 , m 0 , some of the variables will become zero or equal to each other. One can enumerate all possible relations of this kind: Probably there are some other conditions leading to the reduction in the number of variables, but we will restrict ourselves only to the conditions given in equation (3.18). Strictly speaking, at the expense of two parameters one can fulfill two conditions from (3.18). Fulfilling two conditions from the list (3.18) does not ensure simultaneous reduction in the number of variables for all three functions in equation (3.15). However, we will try to find whether it is still possible to choose two arbitrary parameters and fulfill three conditions from the list (3.18). Out of 33 equations from the list (3.18) we created 5456 different systems of equations with 3 equations in each system. Solutions of these systems of equations and analysis of these solutions were performed using MAPLE. It turns out that for some values of m 2 0 , x 3 three conditions from the list (3.18) were fulfilled. In particular, one of such solution reads where The definition of r ijk is given in Appendix A. For x 3 , m 2 0 given in (3.19) the following conditions were fulfilled: Substituting these values of s j0 and m 0 into equation (3.15), we find Explicit expressions for r (i) jkl are given in Appendix A. Thus, the solution of the functional equation (3.15) at m 2 1 = m 2 2 = m 2 3 = 0 is a sum of three terms, each of which is proportional to the same integral ξ (d) 3 depending on two variables different for each term. We notice that µ ijk is a kind of effective mass depending on kinematic invariants.

Verification of the solution of the functional equation
The obtained expression for the integral I (d) 3 given in equation (3.22) is in fact the solution of the equation (3.12). In order to prove this, we substitute the expression for I   with three variables in terms of functions depending on two variables. As we will see in the next subsection, the integral ξ (d) 3 is easier to evaluate than the initial integral. It should be noted that both the Feynman parameter representation for the integral ξ (d) 3 and the dimensional recurrence relation are simpler for this integral than for the initial integral I (d) 3 .
To conclude this section, we consider a particular case of equation (3.22), namely, the case when one of the kinematic invariants, say s 12 , is equal to zero. In this case, is a linear combination of three massless functions ξ (d) The expression for ξ 3 (0; s ij ) can be obtained from the recurrence relation with respect to the spacetime dimension d, which one can find, for example, in Ref. [6]: We draw attention to the coefficient in front of the small imaginary part iη in equation (3.28). The value of this coefficient is important for the analytic continuation of the results in those cases when I  [35], [36].
Setting m 2 = 0 in equation (3.26) and taking into account (3.27), yield Substituting (3.30) in (3.25), we find: Evaluating the Feynman parameter integral for the I (d) 3 gives the same result. Expression (3.31) will be used in the next section for calculating the box type integral.

Analytic evaluation of the integral ξ (d) 3
The analytic expression for the integral ξ  32) or by solving dimensional recurrence relation (3.26) with respect to d. We prefer to use the latter method because it is easier to keep the trace of the term iη at all steps of derivation. To solve the dimensional recurrence relation, we employ the method described in Ref. [7]. For |q 2 | ≤ |4m 2 | the solution of the recurrence relation (3.26) reads where C 3 (q 2 , d) is a periodic function C 3 (q 2 , d) = C 3 (q 2 , d + 2). For this function one can obtain a differential equation with respect to q 2 . It can be derived from the differential equation for ξ 3 (m 2 , q 2 ). In Ref. [37], it was shown that the derivatives with respect to kinematic variables for the Feynman integrals can be written in terms of integrals with shifted space-time dimension d and additional powers of propagators. In our case, the derivative reads: , (3.34) where the values of the kinematic invariants s ij = (p i −p j ) 2 must correspond to those of the integral I (d) 3 in (3.23). With the help of the recurrence relations presented in Ref. [6] the integral on the right-hand side of (3.34) can be reduced to the set of basis integrals. This reduction results in the first-order inhomogeneous differential equation: Note that the term with ξ ln 2m 2 + q 2 + (q 2 + 4m 2 )q 2 depends on a constant of integration K, which may be fixed from the comparison of equation (3.33) with the value of ξ Since ξ (d) 3 (m 2 ; 0) and the term with 2 F 1 at q 2 = 0 are finite, the contribution proportional to the periodic function C 3 (q 2 , d) must also be finite. Taking the limit q 2 → 0 we find that the contribution coming from the term with C 3 (q 2 , d) will be finite if (3.40) Substitution of (3.40), (3.38) in (3.33) yields This expression is valid for | q 2 /4/ m 2 | < 1. The results for ξ for various kinematic regions as well as its ε = (4 − d)/2 expansion was considered in numerous papers [17,[39][40][41][42][43][44][45]. To obtain ε expansion of the integral I (d) 3 depending on three variables, in our approach we need to know ε expansion of the integral ξ (d) 3 depending only on two variables. In order to get the first term in the ε expansion of the integral ξ (d) 3 , the hypergeometric function 2 F 1 in (3.41) must be expanded up to the order O(ε 2 ). Using the HypExp package [20], we find where The leading term in ε for the integral I (d) 3 (0, 0, 0; s 23 , s 13 , s 12 ) can now be obtained by substituting (3.44) into equation (3.22). The resulting expression has been checked for several values of kinematic variables by numerical program [38] for evaluating loop integrals. We found complete agreement for both space-like and time-like values of the kinematic variables s ij .

Functional reduction of the integral I (d) 4
In this section we will consider analytic evaluation of the one-loop integral with massless internal lines associated with the Feynman diagram with four external legs. At d = 4 the analytic result was presented in Refs. [46][47][48]. For particular values of kinematic variables, analytic evaluation of this integral for arbitrary d was considered, for example, in Refs. [35,36,49,50]. Until now, in d-dimensions the analytic result for this integral for all external legs fully off-shell has been known in terms of multiple hypergeometric series [9], [10] evaluated by using the Mellin-Barnes technique. In the present paper, such a result is derived using a combination of functional equations and dimensional recurrence relations.

Derivation of functional equation for the integral I (d) 4 and its solution
In the present paper, we will be concerned with the solution of the functional equation for the box integral with massless internal propagators. However, to find a solution for the integral with all massless internal propagators, we will use the functional equation for the integral with massive internal lines. For this reason, we start our consideration with the functional equation for the box integral with all massive internal lines where The functional equation for this integral may be obtained from the algebraic relation [13]: Here the momentum p 0 is a combination of p 1 , ..., p 4 and m 2 0 , x j must satisfy the following conditions: The solution of this system for x 1 ,x 4 is where Λ 4 is a root of the quadratic equation,   where z 2 , z 3 are arbitrary parameters, q i0 are defined as: q 10 = −λ 4 q 14 + z 3 (q 13 − q 14 ) + z 2 (q 12 − q 14 ) + q 14 , q 20 = (q 12 − q 24 )λ 4 + z 3 (q 23 − q 24 ) − z 2 q 24 + q 24 , and λ 4 is the solution of the quadratic equation: with In order to find a solution of equation (4.12), we will exploit another functional equation, a more general one. Such an equation will be obtained from equation In a manner similar to that for the integral I   where B The explicit expressions for r jkls , r At the next step we will try to reduce integrals B

(d)
ijk with four variables to a combination of integrals with fewer variables. To achieve this goal, we will again use equation (4.10). Setting in this equation Here the parameters x r , m 2 0 are required to obey the following conditions: In order to find conditions on arbitrary parameters for which the number of variables in all integrals on the right-hand side of (4.24) is simultaneously diminishing, we will follow the same strategy which was employed for the functional reduction of the integral I Substituting these values into equation (4.24), we get (4.28) and µ ijk is defined in equation (3.20). The explicit expressions for r (4.29) We have not found relationships reducing ξ (d) 4 to simpler integrals with fewer arguments. Thus, using the two step functional reduction, we expressed the integral depending on six variables in terms of integrals ξ  will be presented in subsection 4.4.

Verification of the solution of the functional equation
Now we will show that the obtained expression (4.29) is the solution of the functional equation (4.12). Substituting I from (4.29) into the right-and left-hand sides of relation (4.12) we obtain 60 terms. The arguments of these functions are the ratios of rather big polynomials containing various powers of square roots of some other polynomials. However, after complicated simplifications of these arguments they became rather compact and it turns out that the situation regarding the integral I integrals, 40 terms with a rather nontrivial dependence on the parameters z 2 , z 3 , cancel each other. The remaining 10 terms were exactly canceled by 10 terms from the left-hand side. is needed for some kinematic variables s ij equal to zero [48], [35,36,[51][52][53]. For this reason it would be interesting to study possible simplifications of relation (4.29) for these particular values of the kinematic variables. In this section, we will use the shorthand   .

Reduction equations for I
(4.33) The considered case is the most complicated one. The number of terms in (4.32) is less than in the general case but the remaining integrals are of the same complexity.
b) The case s 12 = s 34 = 0. This case is simpler than the previous one.     . (4.35) As we will see in the next subsection each integral on the right-hand side of (4.34) is a combination of Gauss hypergeometric functions.
We see that in all the cases but a) the integrals I are combinations of integrals ξ (d) 4 (0, µ 4 , s ij ) with different arguments. The analytic expression for this integral will be given in the next subsection.

Analytic evaluation of the integral ξ (d) 4
An analytic result for the integral ξ where However, we prefer to exploit the recurrence relation with respect to the space-time dimension d. There are several reasons for using this method. First, it is easier to keep the trace of the small imaginary term iη; and second, the resulting expressions usually have a rather compact form in terms of rapidly converging hypergeometric series.
The dimensional recurrence relation for the integral ξ (d) 4 can be obtained from the results for the integral with the general kinematics given in Refs. [6], [54] The solution of this equation can be obtained by using the method described in Ref. [7] and it reads Differentiating relation (4.43) with respect to s ij and taking into account equation (3.35), we obtain On the other hand, we can write this derivative of ξ (d) 4 with respect to s ij in terms of the d + 2 dimensional integral with additional powers of propagators [37] ∂ ∂s ij ξ (d) where P 1 = (k 1 − p 1 ) 2 + iη, P 2 = (k 1 − p 2 ) 2 + iη, The kinematic invariants s ij = (p i − p j ) 2 in (4.46) are to be the same as those for the integral (4.28). After applying the recurrence relations [6], [54] to reduce this integral to a set of basis integrals, we obtain for ξ (d) 4 the first-order differential equation Using this relation, we find with equal masses is well known [55], [56] A combination of equations (4.44), (4.47) yields the first order differential equation for C 4 (s ij , d): The solution of equation (4.52) is where κ c is a constant of integration, and The constant of integration κ c is easy to fix from the finiteness condition for C 4 (s ij , d) as s ij → 0: Finally, combining all the contributions, we find The analytic result for the first sum in equation (4.43) can be obtained with the help of equation (3.41), and it is (4.60) Substituting (4.60), (4.59) into equation (4.43) yields It is important to note that for evaluating analytic results for all the three particular cases the only function, namely ξ (d) 4 (0, µ 4 ; s ij ), is required. This is one of the advantages of our functional reduction method. for particular values of variables In this subsection, we will present the leading term in the expansion of the integral I for general kinematics for a future publication.
To obtain ε = (4 − d)/2 expansion of I (d) 4 for the cases b), c) and d), we need to know the first terms in the expansion only for the integral ξ (d) 4 (0, µ 4 , s ij ). Plugging expansions for the 2 F 1 function, obtained by using the HypExp package [20], and for the integral ξ where L ij = ln(− s ij ). (4.68) We compared numerical values derived from this formula with the results obtained by the numerical program SecDec [38] and found complete agreement for both the space-and time-like values of the scalar invariants s ij . Using equation (4.67), it is easy to derive from equation (4.38) the leading term in the ε expansion for the integral I  Here and in formulae below we used L ij defined in equation (4.68).
For the case c) the leading term in ε is:  Integrals appearing after the last iteration of the functional reduction were evaluated by using the dimensional recurrence relations. A distinctive feature of these recurrence relations is that they are the first order inhomogeneous equations and the inhomogeneous part has only one term. This significantly simplified their solution.
There are many directions for future applications of the proposed method. First of all, we are going to apply our approach to the reduction of massless one-loop scalar integrals associated with diagrams with five and six external legs. Also, the method can be extended without problems to the one-loop integrals with massive internal propagators.
Another important direction for future research will be the extension of the method of functional reduction to the evaluation of multiloop integrals. One can easily write down functional equations for multiloop integrals, but to elaborate a systematic algorithm one should solve a number of problems. At the present time, we are working on the solution of these problems.

Appendix A
In this appendix, we give the definition of the Gram determinants and explicit formulae for polynomials occurring in equations (3.22), (4.20), (4.29).
where s 2 ij = (p i − p j ) 2 , p i are combinations of external momenta flowing through i-th lines, respectively, and m i is the mass of the i-th line. Where no confusion can arise, we simply refer to the above functions as ∆ n , G n−1 . We will also use an indexed notation for ∆ n and G n−1 . . , p in ). (7.5) Our results depend on the ratios of λ i 1 i 2 ...in and g i 1 i 2 ...in and, therefore, it is convenient to introduce the notation The imaginary part of r is rather simple. Using one shows that to all orders in η ..in η, (7.8) and, therefore, the causal η prescription for r is (with the same η for all masses) r ij...k | m 2 j −iη = r ij...k | m 2 j − iη. (7.9) For the reader's convenience we present below explicit expressions for the ratios of Gram determinants and their derivatives occurring in the reduction formulae for the integrals I and I (d) 4 : [2s ij s ik s kn − 2s ij s 2 kn + 2s ij s kn s in − 2s 2 ik s jn + 2s ik s jk s in + 2s ik s jn s kn + 2s ik s jn s in − 4s ik s kn s in + 2s jk s kn s in − 2s jk s 2 in ], r (k) ijkn = 1 g ijkn [2s ij s ik s jn − 2s 2 ij s kn + 2s ij s jk s in + 2s ij s jn s kn − 4s ij s jn s in + 2s ij s kn s in − 2s ik s 2 jn + 2s ik s jn s in + 2s jk s jn s in − 2s jk s 2 in ], [2s ij s ik s jn − 2s 2 ij s kn − 4s ij s ik s jk + 2s ij s ik s kn + 2s ij s jk s kn + 2s ij s jk s in − 2s 2 ik s jn + 2s ik s jk s jn + 2s ik s jk s in − 2s 2 jk s in ], r ljkn = − λ ljkn g ljkn = r ljkn − iη, (7.11) where g njk = 2(s 2 nj + s 2 nk + s 2 jk − 2s nj s nk − 2s nj s jk − 2s nk s jk ), λ njk = −2s nj s jk s nk + iηg njk , (7.12) g ljkn = 4s 2 lj s kn + 4s lj s lk s jk − 4s lj s lk s jn − 4s lj s lk s kn − 4s lj s ln s jk + 4s lj s ln s jn − 4s lj s ln s kn − 4s lj s jk s kn − 4s lj s jn s kn + 4s lj s 2 kn + 4s 2 lk s jn − 4s lk s ln s jk − 4s lk s ln s jn + 4s lk s ln s kn − 4s lk s jk s jn + 4s lk s 2 jn − 4s lk s jn s kn + 4s 2 ln s jk + 4s ln s 2 jk − 4s ln s jk s jn − 4s ln s jk s kn + 4s jk s jn s kn , (7.13) λ ljkn = −2s lj s lk s jn s kn − 2s lj s ln s jk s kn − 2s lk s ln s jk s jn + s 2 lj s 2 kn + s 2 lk s 2 jn + s 2 ln s 2 jk + iηg ljkn . (7.14) 8 Appendix B In this appendix, we describe derivation of the analytic dependence of the one-loop massless propagator integral (3.27) on the small imaginary part iη added to the propagators. Expression This expression agrees with the one presented in Ref. [56]. For the massless case in expression (8.3) we must set m 2 = 0 or As η → 0, the second term in equation (8.3) is exponentially small compared to the first term and, therefore, it can be neglected. Thus, the leading contribution to the integral I We simplified the first term in (8.3) by dropping terms proportional to η, keeping the structure of the branch point in the vicinity of q 2 = 0. The leading term in (8.5) is in agreement with (3.27).

Appendix C
For the sake of completeness, we present in this appendix useful formulae for the Appell F 1 and Gauss 2 F 1 hypergeometric functions.