{'fraction_non_alphanumeric': 0.04871683340582862, 'fraction_numerical': 0.024081616513108466, 'mean_word_length': 4.021842732327244, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The massless two-loop two-point function

We consider the massless two-loop two-point function with arbitrary powers of the propagators and derive a representation, from which we can obtain the Laurent expansion to any desired order in the dimensional regularization parameter eps. As a side product, we show that in the Laurent expansion of the two-loop integral only rational numbers and multiple zeta values occur. Our method of calculation obtains the two-loop integral as a convolution product of two primitive one-loop integrals. We comment on the generalization of this product structure to higher loop integrals.


Introduction
Quantitative predictions from perturbation theory are crucially linked with our ability to calculate loop integrals. An object of extensive study has been the master two-loop two-point function, with massless internal propagators but arbitrary powers of the propagators. The name "the master two-loop two-point function" is justified, since all other two-loop two-point topologies can be obtained from this one. Allowing arbitrary powers of the propagators is important for threeor four-loop calculations, where the integration over self-energy insertions on the propagators shifts the power of the propagator from unity to 1 + ε. (As usual we work within dimensional regularization with D = 4 − 2ε.) In this paper we consider the massless two-loop two-point function. This integral is not only of practical importance from a phenomenological perspective, but received also quite some interest from number theorists. As far as the phenomenological side is concerned, this integral is implemented into the MINCER package [1,2] and used for example for the calculation of the total hadronic cross-section in electron-positron annihilation. From the number theoretical perspective there has always been the open question which types of (transcendental) numbers appear in the ε-expansion of this integral. From explicit calculations it is known, that in the lowest orders multiple zeta values occur.
The history for the calculation of this integral dates back, to the best of our knowledge, to 1980, when the ε 0 -coefficient for a specific (non-trivial) combination of powers of propagators was calculated using the Gegenbauer polynomial x-space technique [3]. In 1984 Kazakov obtained for all powers of the propagators of the form ν j = 1 + a j ε the result up to ε 3 , a year later the ε 4 -term followed [4,5,6]. The momentum of pushing the ε-expansion of this integral further was then taken up by Broadhurst and collaborators: In 1986 the ε 5 -term was calculated, in 1988 followed the ε 6 term [7,8]. In the mean time the original Gegenbauer technique had been refined [9] and a representation of the integral for a specific combination of the powers of the propagators in terms of hypergeometric 3 F 2 -functions with unit arguments was found [9,10]. The year 1996 brought an advance of two further terms in the expansion (ε 7 and ε 8 ) [10]. Finally, last year the ε 9 -term was announced [11]. The calculation of most of these terms relied on symmetry properties of the two-loop integral. It is known that this two-loop integral has the symmetry group Z 2 × S 6 [12,7,8], which is of order 1440, e.g. there are 1440 symmetry relations. The calculations cited above exploited the fact that the symmetry properties allow to reconstruct the result of the expansion up to order ε 9 for powers of the propagators of the form ν j = 1 + a j ε from the result of the expansion, where two adjacent propagators occur with unit power. However, it is known that beyond this order the symmetry relations are not sufficient to determine the solution [10]. In view of the last point, the calculation of this two-loop integral cannot be considered to be solved to a satisfactory level. We quote from the latest publication of Broadhurst [11]: "It is one of the many scandals of our limited understanding of the analytical content of perturbative quantum field theory that, despite many years of intense effort, we still do not know whether multiple zeta values suffice for even the Taylor expansion of the two-loop integral." In this paper we proof that multiple zeta values are sufficient (theorem 4.1). This is one of the main results of this paper and solves a long standing open problem. We also show how to obtain the ε-expansion to arbitrary order for arbitrary powers of the propagators by calculating this integral with the help of a new method. For an introduction to calculational techniques for multi-loop integrals we refer to [13,14,15].
The second important result of this paper is the method we employed for the calculation of the two-loop integral: We obtain the two-loop integral as a convolution product of two primitive one-loop integrals. The convolution product can be evaluated in terms of nested sums [16,17], and the ε-expansion of the original integral is obtained from the ε-expansion of these sums. Generalizations of this product structure can be useful for the extension of exisiting packages for two-point functions (like MINCER [1,2]) to four loops. In this paper we only briefly comment on this possibility and focus on the two-loop two-point function. In view of later applications for the calculations of higher orders we show that our method for the calculation can be implemented efficiently on a computer. This paper is organized as follows: In sect. 2 we define the two-loop integral and summarize known facts about this integral. In sect. 3 we show that the two-loop integral can be written as a convolution product of two one-loop integrals. We also discuss the factorization for three-loop two-point functions. In sect. 4 we use this product structure to obtain the ε-expansion for the two-loop integral. In sect. 5 we report on the implementation of our formulae into a symbolic computer code. Sect. 6 contains a summary and our conclusions. An appendix collects some useful formulae for integral transformations.

Review of known results for the two-loop integral
The object of investigation is the following five-propagator integral I (2,5) The superscripts for an integralÎ (l,n) indicate the number of loops l and the number of propagators n. The prefactor in front of the integral is inserted for later convenience. It ensures thatÎ (2,5) is independent of (−p 2 ) and avoids a proliferation of Euler's constant γ E . The corresponding Feynman diagram for this two-loop integral is shown on the l.h.s of fig. 1. It is for our purpose sufficient to assume that the exponents ν j are of the form where the n j are positive integers and the a j are non-negative real numbers. The ε-dependence in the powers of the propagators arises from higher loop integrals by integrating out simple one-loop two-point insertions into the propagators of eq. (1).
Integration-by-part identities [18,19,20] relate integrals with different powers of the propagators. For example, from the triangle rule [21] we obtain the identities Here, the operators i + and i − raise, respectively lower, the power of propagator i by one. The integration-by-part relations can be used to relate the integral where all propagators occur to a positive integer power to simpler topologies, e.g. where one propagator is eliminated. The integral in eq.(1) has the obvious symmetries However, there are more symmetries, which relate the integral to itself, up to prefactors of products of Gamma functions. To discuss the symmetry properties of eq. (1) it is convenient to introduce the function F(ν 0 , ν 1 , ν 2 , ν 3 , ν 4 , ν 5 ) related toÎ (2,5) bŷ where and Here and in the following we will use the short-hand notation like ν i jk = ν i + ν j + ν k to denote sums of indices if {i, j, k} ∈ {1, 2, 3, 4, 5}. The function F(ν 0 , ν 1 , ν 2 , ν 3 , ν 4 , ν 5 ) is invariant under the symmetry group Z 2 × S 6 [8]. The symmetric group S 6 is generated by the six-cycle and the transposition The group Z 2 is generated by the reflection Note that in general the generators of the symmetry group do not conserve the positivity a j ≥ 0 in the parameterization of eq. (3).

The product structure of the two-loop integral
In this section we show that the two-loop integral in eq. (1) can be written as a (convolution) product of two one-loop integrals. To this aim we define the following two one-loop integrals: . Note that the integral I (1,3) depends on the kinematic variables −p 2 , −k 2 1 and −k 2 4 only through the dimensionless ratios x and y. The integralÎ (1,2) is easily computed aŝ The one-loop triangle integral I (1,3) can be written as a double Mellin-Barnes representation: where the functionÎ (1,3) (m − ε, ν 2 , ν 3 , ν 5 ; τ, σ) is given bŷ The integration contours are parallel to the imaginary axis, with indentations, if necessary, to separate the "UV"-poles ( x, y): From eq. (13) and eq. (14) one obtains the two-loop integral aŝ In eq. (17) the two-loop integral is obtained as a (double) convolution product of two one-loop integrals. This is pictorially shown in Fig. (1). At the level of Feynman diagrams the product structure is given as the insertion of a Feynman diagram into another Feynman diagram. It is a well known fact that convolution products can be turned into ordinary products by applying a suitable integral transformation. In the case at hand, eq. (17) factorizes by performing two inverse Mellin transformations in ν 1 and ν 4 . If one sets one obtains the double inverse Mellin transform of the two-loop integralÎ (2,5) as the product of the one-loop integral I (1,3) with the double inverse Mellin transform of the one-loop integral I (1,2) : Eq. (17) or eq. (19) is the advertised factorization of the two-loop integral into two one-loop integrals. At the level of Feynman diagrams the product is similar to the insertion operation defined by Kreimer within the context of renormalization [22,23,24]. In the context discussed here, insertions occur only at one specified place. This implies that the product is associative. Note that in general an insertion product, which allows insertions at several places, is not associative. The factorization property, e.g. that non-primitive graphs can be written as convolution products of primitive graphs, generalizes to higher loops. For the three-loop two-point functions there are three basic topologies, usually named the ladder ("LA") topology, the Benz ("BE") topology and the non-planar ("NO") topology. They are shown on the l.h.s of fig. (2). The ladder and Benz topologies are given as convolution products of three one-loop graphs: Figure 2: Factorization at three loops: The first two topologies factorize into one-loop diagrams, whereas the last topology factorizes into a one-loop diagram and a two-loop diagram. A graph on the right side of the product operator ⋆ is inserted into the shaded vertex of the graph to the left of the product operator.
The momenta are defined for the ladder topology by For the Benz topology we havê The momenta are defined for the Benz topology by The non-planar topology factorizes only into a one-loop two-point function and the crossed two-loop three-point function: Here,Î of the two-loop crossed vertex function The momenta for eq. (22) and eq. (24) are defined by k 4 = k 3 − p, k 5 The variables x and y are defined by x = (−p 2 )/(−k 2 1 ) and y = (−p 2 )/(−k 2 6 ). For all topologies the factorization is shown pictorially in fig. (2).

Evaluation of the two-loop integral
Eq. (17) is the starting point for the further evaluation of the two-loop integral: The strategy for the evaluation is as follows: We first close the contours and evaluate the integrals with the help of the residuum theorem. This technique has a long history, a recent example is the calculation of the triple box [25]. The semi-circles at infinity needed to close the contours give a vanishing contribution provided that From the residues we immediately obtain nested sums. All these sums can be brought to a standard form and then expanded into a Laurent series in ε with the help of the algorithms described in [16]. We use the program "nestedsums" [17] and the "GiNaC"-library [26] for this purpose.
The algorithms A and B of [16] are generalizations of algorithms for the manipulations of harmonic sums and harmonic polylogarithms, described in [27,28,29]. Closing the contour of the σ-integration to the right, one picks up the residues of Γ(−σ), Γ(−σ + m − ε − ν 35 ) and Γ(−σ − τ − m + ε + ν 14 ). All other σ-dependent Gamma functions in the numerator have poles to the left of the contour and therefore do not contribute. The basic formula for the residuum of Euler's Gamma-function reads: The two-loop integral can therefore be written as a sum of three terms, where each term is obtained by taking the residues of one Gamma function from the set Γ(−σ), . Explicitly, with , The procedure is then repeated by closing the τ-integration contour to the right. For example, for the term T (1) we have to evaluate the residues of Γ(−τ), Γ(−τ + m − ε − ν 25 ) and Γ(−τ − n − m + ε + ν 14 ). The residues of Γ(τ + m − ε − ν 1 ), Γ(τ + n − m + ε + ν 235 ) or Γ(τ + n + ν 5 ) are always to the left of the contour. In summary, the two-loop integral can be written as the sum of the following terms: Table (1) shows the correspondence between each term T (k,l) and the Gamma functions, from which the residues are taken. The term T (3,3) corresponding to the residues in Γ(−τ + n + ν 1 ) gives no contribution. This Gamma function is always accompanied by 1/Γ(−τ + ν 1 ). Since is free of poles we have Table 1: Correspondence between the terms T (k,l) and the Gamma functions, from which the residues are taken.
Furthermore, we can show that the terms T (1,3) and T (3,4) as well as the terms T (2,3) and T (3,5) cancel each other: To present the results after all residues have been taken, we introduce two functions G ± with ten arguments each: together with two operators L d and R d acting on the arguments as follows: ThenÎ (2,5) Here, 1 denotes the identity operator with a trivial action on the arguments of the functions G ± .
To proceed further, we now show how to transform the functions G ± to a standard form, such that they can be expanded in ε with algorithms of [16,17]. Since the shift operators L and R modify only the arguments but not the structure of the functions G ± , it is sufficient to discuss one example for each function. For the function G + we discuss as an example the first term (without any application of the shift operators L or R ) of eq. (38), which corresponds to the term T (1,1) .
If ν 1 and ν 4 are not positive integers, one may use the reflection formula to flip the Gamma functions, where a summation index occurs with a negative sign. One obtains in that case .
This is a double infinite sum with unit arguments and can be evaluated with algorithm B of [16,17]. The sum in eq.(41) is a generalization of the first Appell function. The Laurent expansion to any fixed order will contain only rational numbers and multiple zeta values. If ν 1 or ν 4 are positive integers, the infinite sums over j or n terminate. For example, if ν 1 is a positive integer, while ν 4 is not, one obtains .
This is a finite sum of single infinite sums, more precisely, it is a sum of ν 1 terms, each containing a 4 F 3 hypergeometric function with unit argument. These are evaluated with algorithm A of [16,17]. Again, the Laurent expansion to any fixed order will contain only rational numbers and multiple zeta values. The case where ν 4 is a positive integer while ν 1 is not, is completely analog. If both ν 1 and ν 4 are positive integers, one obtains e.g. a finite sum of products of Gamma functions. This is rather trivial and the Laurent expansion contains only rational numbers and zeta values. We now turn to the function G − . As an example we discuss If ν 1 is not a positive integer, we split the summation region into the regions j ≤ n and n ≤ j and subtract the double counted diagonal j = n as in We obtain The sums are of the same type as discussed in eq. (41) and eq. (42), therefore the Laurent expansion contains only rational numbers and multiple zeta values. If ν 1 is a positive integer, we observe that in eq. (44) we have the building block Therefore we obtain a non-vanishing contribution only if this sequence does not contain zero. This is the case for Region I: j < ν 1 , Region II: j > n + ν 1 − 1.
Region I contains only a finite number of residues from Γ(−τ), whereas region II contains an infinite sum. We obtain Once again, the sums are of the same type as discussed in eq. (41) and eq. (42), and the Laurent expansion contains only rational numbers and multiple zeta values.
To summarize, we are able to write the two-loop integralÎ (2,5) as a combination of terms, which can be expanded in ε to arbitrary order with the help of algorithms A and B of [16,17]. All these terms occur with unit arguments, therefore the Laurent expansion ofÎ (2,5) involves only rational numbers and multiple zeta values. In particular, in the representation derived here, no alternating Euler sums occur. This answers a question raised recently by Broadhurst [11]: Theorem 4.1 Multiple zeta values are sufficient for the Laurent expansion of the two-loop inte-gralÎ (2,5) (m − ε, ν 1 , ν 2 , ν 3 , ν 4 , ν 5 ), if all powers of the propagators are of the form ν j = n j + a j ε, where the n j are positive integers and the a j are non-negative real numbers.
We have shown above, that the Laurent expansion of the two-loop integralÎ (2,5) can be expressed in multiple zeta values if the inequalities eq. (27) are satisfied. These inequalities ensure that the semi-circles at infinity give a vanishing contribution. If for a specific combination of powers of propagators these inequalities are not satisfied, one may use the integration-by-part identities eq. (4) and the symmetry relations eq. (5) and eq. (9) - (11) to express the original integral as a linear combination of integrals, which fullfill the inequalities. The coefficients of this linear combination are either rational numbers or Gamma functions, which in turn expand into zeta values. Therefore the theorem holds also in the general case.
The restriction to a j ≥ 0 in the parameterization ν j = n j + a j ε ensures that the different cases we treated above (e.g. whether the argument of a Gamma function is an integer or not) are sufficient. This condition could be relaxed at the expense of a more extensive case study. However, in all practical calculations the integrals occur with a j ≥ 0.
Finally we remark on a technical detail: Our method of calculation expresses the two-loop integral as a combination of several sums with unit arguments. Although the final result is finite, it is not guaranteed that all individual sums are convergent. Furthermore, the algorithms used for the Laurent expansion rely on partial fractioning. This may split a convergent sum into two divergent pieces, as illustrated by the example of the convergent sum Partial fractioning splits 1 n(n + 1) The problem is easily circumvented by introducing a finite upper summation limit N or by multiplying the summand by x n and taking the limit x → 1 in the end. With both methods, the divergent pieces cancel at the end of the day and one obtains the correct and finite result. Within the second method, one inserts into the integrand of eq. (25). Then all sums are convergent, provided 0 ≤ x 1 < 1, 0 ≤ x 2 < 1 and x 2 < x 1 . The Laurent expansion for this more general expression will contain multiple polylogarithms in x 1 and x 2 . The result for the two-loop integralÎ (2,5) is recovered by first taking the limit x 1 → 1 and then the limit x 2 → 1. Note that the order of the limits cannot be exchanged. It is easily seen that in the limits the multiple polylogarithms reduce to multiple zeta values. 0 1 2  3  4  5  6  7  weight  3 4 5  6  7  8  9  10  time  29 54 99 185 375 910 2997 11741  memory 6 7 8 12  30 104 397 1970   Table 2: CPU time in seconds and required memory in MB for the expansion ofÎ (2,5) (2 − ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε) up to the indicated order/weight on a PC (1.6 GHz Athlon with 2 GB RAM).

Results, checks and performance
We have implemented the results of the previous section into a C++ program, which calculates symbolically the Laurent expansion of the two-loop integral for a user-specified set of parameters (m, ν 1 , ..., ν 5 ) up to the desired order in ε. The program uses the "nestedsums"-library [17] and the "GiNaC"-library [26]. To simplify our results we use the Gröbner basis for multiple zeta values provided in [30].
To check our implementation we have written two independent programs and verified that they agree. Further we have checked that for trivial cases (e.g. all powers of the propagators integers) we obtain the correct (known) result. In addition, we have verified the symmetry relations eq.(5). Finally we compared for non-integer powers of the propagators with known results from the literature. As an example we quote the result for one particular integral: (1 − 2ε)Î (2,5) (2 − ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε) = 6ζ 3 + 9ζ 4 ε + 372ζ 5 ε 2 + 915ζ 6 − 864ζ 2 3 ε 3 + (18450ζ 7 − 2592ζ 4 ζ 3 ) ε 4 + 50259ζ 8 − 76680ζ 5 ζ 3 − 2592ζ 6,2 ε 5 + 905368ζ 9 − 200340ζ 6 ζ 3 − 130572ζ 5 ζ 4 + 66384ζ 3 3 ε 6 + 2955330ζ 10 − 68688ζ 8,2 − 3659904ζ 7 ζ 3 − 1777680ζ 2 5 + 298728ζ 4 ζ 2 For multiple zeta values we use the notation We have compared our result up to weight 10 with [10] and found agreement. The order to which we can expand the two-loop integral in the parameter ε is not limited by our method of calculation. The only restriction arises from hardware constraints (e.g. available memory and CPU time). Table (2) shows for the case of the integralÎ (2,5) (2 − ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε, 1 + ε) the dependence of the required CPU time and the required memory on the order of ε to which the Laurent expansion is calculated. In table (2) we also indicate the highest weight of multiple zeta values, which occur within a given order. Note that the weight is a more accurate measure of the complexity of a calculation, since in the example discussed here individual terms start at 1/ε 3 , but in the sum the coefficients of the terms 1/ε 3 , 1/ε 2 and 1/ε vanish. In general we expect terms up to weight 2l − 1 to occur in the finite part of a l-loop two-point function. We see from table (2) memory needed to calculate each of these integrals is slightly below the corresponding numbers indicated in table (2).

Summary and conclusions
In this paper we calculated the ε-expansion of the massless master two-loop two-point function with arbitrary powers of the propagators. We showed that to all orders in ε rational numbers and multiple zeta values are sufficient to express the result. Our method of calculation obtained the two-loop integral from a convolution of two one-loop integrals. We also discussed the corresponding factorization for three-loop two-point functions. Finally we demonstrated that our method can be implemented efficiently on a computer.
Then the Mellin transform is defined for c 0 < Re σ < c 1 by The inverse Mellin transform is given by The integration contour is parallel to the imaginary axis and c 0 < Re γ < c 1 . There is a close relation between the Laplace and the Mellin transform: If f (t) = h(e −t ), then f L (σ) = h M (σ). Note that the inversion formulas can be obtained from the known inversion formula for the Fourier transform.
To a function f (t) we associate the Fourier transform The inverse transform is given by The Laplace transform and the Fourier transform of a function f (t) are related by As an example for the Mellin transform we consider the function with Mellin transform h M (σ) = Γ(−σ)Γ(σ + c)/Γ(c). For Re(−c) < Re γ < 0 we have We often deal with integrals of the form I = 1 2πi γ+i∞ γ−i∞ dσ Γ(σ + a 1 )...Γ(σ + a m ) Γ(σ + c 2 )...Γ(σ + c p ) If max (Re(−a 1 ), ..., Re(−a m )) < min (Re(b 1 ), ..., Re(b n )) the contour can be chosen as a straight line parallel to the imaginary axis with max (Re(−a 1 ), ..., Re(−a m )) < Reγ < min (Re(b 1 ), ..., Re(b n )) , otherwise the contour is indented, such that the residues of Γ(σ + a 1 ), ..., Γ(σ + a m ) are to the right of the contour, whereas the residues of Γ(−σ + b 1 ), ..., Γ(−σ + b n ) are to the left of the contour. We further set Then the integral eq. (68) converges absolutely for α > 0 [33] and defines an analytic function in |arg x| < min π, α π 2 . (71) The integral eq. (68) is most conveniently evaluated with the help of the residuum theorem by closing the contour to the left or to the right. Therefore we need to know under which conditions the semi-circle at inifinty used to close the contour gives a vanishing contribution. This is obviously the case for |x| < 1 if we close the contour to the left, and for |x| > 1, if we close the contour to the right. The case |x| = 1 deserves some special attention. One can show that in the case β = 0 the semi-circle gives a vanishing contribution, provided λ < −1.
To derive this result, Barnes asymptotic expansion of the Gamma function for large x is useful: where the Bernoulli polynomials B n (x) are given in terms of the Bernoulli numbers B j by B n (x) = n ∑ j=0 n j B j x n− j .