{'fraction_non_alphanumeric': 0.03530858838000145, 'fraction_numerical': 0.011142464365820129, 'mean_word_length': 4.2763886237831645, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of Massless Feynman Integrals using Harmonic Sums

A method for the evaluation of the epsilon expansion of multi-loop massless Feynman integrals is introduced. This method is based on the Gegenbauer polynomial technique and the expansion of the Gamma function in terms of harmonic sums. Algorithms for the evaluation of nested and harmonic sums are used to reduce the expressions to get analytical or numerical results for the expansion coefficients. Methods to increase the precision of numerical results are discussed.


Introduction
One of the most important aims of particle physics is the test of models, particularly the standard model. Testing the standard model means, on one hand, to increase the accuracy of its parameters and, on the other hand, the search for deviations of its predictions from measured quantities to discover new effects. To do this experiments are accomplished with increasing accuracy of measurement. Therefore theoretical calculations have to be performed with increasing numerical precision.
A parameter of the standard model of particular interest is the strong coupling constant α s . It can, for example, be determined via the so-called R-ratio, the normalized total cross section of the process e + e − → hadrons, R(s) = σ(e + e − → hadrons) σ(e + e − → µ + µ − ) .
In order to achieve high precision one has, within perturbative quantum field theory, to calculate high orders of the perturbation series. This means, we need to calculate large numbers of more and more complicated Feynman diagrams. Therefore different methods for the evaluation of such multiloop Feynman diagrams have been developed. (See e.g. [1,2,3].) The quantity R(s) is nowadays known to O(α 3 s ) precision [4,5] and first results in O(α 4 s ) have been published [6]. For a review see Refs. [7,8].
A prescription to solve problems with a large number of complicated Feynman diagrams is to reduce them to a set of master integrals, e.g. by means of the traditional integration by parts method, and then evaluate these.
The aim of the project described in this article is to develop a method for the calculation of massless integrals. More precisely, we want to consider the planar (P ) and non-planar (N) three-loop propagator diagrams: They can be considered as the master integrals of the corresponding class of Feynman diagrams, which are "building blocks" of the calculations for R(s) in O(α 4 s ).
We are interested in the ε-expansion of the diagrams, where ε is the parameter of dimensional regularization. This means, the integrals are formally calculated in D = 4 − 2ε dimensions and then expanded in ε. By means of this procedure UV divergencies of diagrams manifest themselves as poles in ε and can then be cancelled by means of renormalization.
It is thus essential for renormalization to know the pole structure of a Feynman diagram and the finite part of the expansion. However, when one has to deal with products of graphs with poles as in the case of the following diagram, which appears within the mentioned calculations of R(s), one has to calculate higher ε orders of the second factor in order to ensure, that the result is correct up to the finite term. This is why we are interested in higher orders of the ε expansion of the diagrams under consideration.
The O(ε 1 ) terms of the diagrams N and P are necessary for the computation of the β function of the scalar φ 4 model at five-loop level. They have been known since long from Ref. [9] and [10].
The calculation of the planar diagram is somewhat simpler than that of the nonplanar one, because it could be related to the generic massless two-loop two-point function by means of integration by parts [11]. The latter was studied in many works. e.g. [10,12,13,14,15,16] and references therein. Recently, in Ref. [17] it was in a sense computed analytically. Thus one can say that the planar diagram is known to arbitrary order and certainly to sufficient precision for practical needs.
These numbers have been used to check the method and to examine the properties of the numerical algorithms described later.
The determination of the ε 2 coefficient of N is one of the aims of the project presented in this paper. An independent evaluation can be found in Ref. [19], the result was given in Ref. [20]. A numerical approach to the calculation of integrals of this type is given e.g. in Ref. [21].
The paper is organized as follows: Section 2 shows how a representation of the diagrams N and P in terms of sums can be found by means of the Gegenbauer polynomial method. Section 3 gives an overview over the Harmonic Sums and algorithms to deal with them. It is shown how the summation algorithms are applied to evaluate the sums representing the diagrams. For the N diagram an analytical solution by means of the described algorithms is not yet possible. In this case one can get numerical results by summing up the series. Section 4 reviews methods to accelerate the convergence of such processes to get better results. By means of these numerical methods it is possible to calculate the expansion coefficients with high precision. For the non-planar diagram this is described in Section 5. The planar topology can be calculated analytically with the help of the summation algorithms, which is demonstrated in Section 6.

The Gegenbauer Polynomial x-Space Technique
The Gegenbauer Polynomial x-Space Technique (GPXT) to evaluate Feynman Diagrams in coordinate space was introduced in Ref. [12]. Applications can be found e.g. in Refs. [22,23], an extension of the technique in Ref. [24]. In this Section we will apply the method to the diagram N. The results for the planar diagram can be obtained analogously. Labeling the momenta as shown in Fig. 1 the integral representation of N reads By means of Fourier transformation where λ = 1 − ε = D 2 − 1, the integrand is transformed to position space. The integrals over the momenta can then be evaluated to delta functions: After the evaluation of the delta functions and choosing an appropriate coordinate system (see Fig. 1, where the root vertex is denoted by 0) the denominators get a simple and uniform shape: Spherical coordinates are introduced by means of where Now the propagator terms and the exponential function can be expanded in Gegenbauer polynomials using formulae (A.8) and (A.11) of appendix A. This leads to a representation of N as a multiple sum The angular integrations can be evaluated using the properties of the Gegenbauer polynomials. Each application of the orthogonality relation (A.4) reduces two C λ n to one and produces a Kronecker delta. Finally, one arrives at the expression The complicated structure of the triple sum is introduced by Eq. (A.14), which is applied when there are two polynomials with the same arguments.
To handle the radial integrals the region of integration is distributed in 120 regions according to the relative size of r 1 , r 2 , r 3 , r 4 , r 5 . In each region the expressions M n λ , Eq. (12), evaluate to rational expressions. For simplification k 2 is set to 1; the k-dependence can later be restored by dimensional arguments. Thus we are lead to integrals like Actually, due to the symmetry of the diagram, only 30 of these 120 terms are different. The calculation of these can be automated. * After performing all integrations and rewriting λ = 1 − ε the integral N is expressed in terms of a triple sum: where R(ε, l, m, n) = 4 (l + 1) (l + m + n + 4) (m − 2ε + 1) (l + m + n − 6ε + 2) + 8 (l + 1) (l + m + n + 4) (l + m + n − 4ε + 2) (l + m + n − 6ε + 2) + . . . (118 similar terms).
Finally, according to Ref. [12], this result can be cast into the so-called G-Form. Rewriting one obtains N S can now be expanded in powers of ε. For each order in ε a triple sum has to be evaluated. The calculation of these sums is the aim of the rest of the paper.

Summation Algorithms
In this Section the Harmonic Sums are introduced and their relationship to the problem under consideration is shown.

Harmonic Sums
Harmonic Sums are defined as [25,26,27] which, for n = ∞, leads to Riemann's zeta function: Higher Harmonic Sums are defined [25] recursively via the relation The number of indices of a harmonic sum is called its depth, the sum of the absolute values of the indices its weight.
A negative index denotes an alternating sum, † In Ref. [25] different algorithms are described to deal with expressions of the form and sums over such. These are implemented in the program package SUMMER written in FORM [28] which is available under the URL http://www.nikhef.nl/~form/FORMdistribution/packages/summer.
The harmonic sums are related to the psi (digamma) function, which is the logarithmic derivative of the Gamma function, and its derivatives, the polygamma functions ψ (n) (x), † Note that this is the convention used in [25]. Other authors denote a factor (−1) m by overlining the index. In this case a negative sign of the index indicates -consequently -a positive power of the corresponding variable.
via the following relations [29]: Using these relations the Gamma function can be expanded in terms of harmonic sums (cf. Ref. [30]) The expansion of the exponential function leads to the expansion formulae used in the following:

Simplification of the sums
The structure of the sum given in Eq. (15), particularly the boundaries of the innermost sum, is very difficult to handle with symbolic manipulation programs. Therefore the first step is the simplification of this sum.
For a given n = ν we have to sum over pairs of indices (l, m), which satisfy Neglecting the condition l + m + n = 2g at first, this is equivalent to which is described by the rectangular area in the (l, m)-plane shown in Fig. 2. Fig. 2. Area of summation for fixed n = ν.
Taking into account the condition l + m + n = 2g this means, that the sum ranges over all combinations of indices in this marked area whose sum is even.
Let us in a next step consider pairs (l, m) with constant sum a = l + m, i.e. points on straight lines parallel to the second bisector. The sum starts with the line m = ν − l (a = ν) and reaches over the shaded rectangle in Fig. 2 to infinity. Taking into account the symmetry of the terms it is sufficient to restrict oneself to the area below or over the first bisector. The contributing points on such a line can be described by the condition k ≤ 2l ≤ k + 2ν, where ν + k = a (cf. Fig. 2).
If we parametrize k = 2κ and k = 2κ + 1, the addends read F (l, ν + 2κ − l, ν) and F (l, ν + 2κ + 1 − l, ν), respectively, and it is obvious, that in the second case the sum of the arguments can never be even, while the terms of the first kind automatically fulfill the condition.
Therefore the original sum is transformed to a much simpler form: If the expansions (29,30) are applied to the sum in Eq. (11), it is represented as a triple sum over rational expressions and products of such and harmonic sums. In the following some algorithms are discussed to simplify such objects.

Rational Expressions
To sum expressions of the form where P (x) is a polynomial in x, one first performs partional fractioning. This leads to simpler summands in which x appears only once in the denominator. The resulting sums can be evaluated easily in terms of zeta or psi functions by means of the formula [31] β x=α If the upper limit is infinity, the second psi function is zero for n > 1.

Algorithms for Nested Sums
In Ref. [26] different algorithms have been published to evaluate nested sums using the properties of higher harmonic sums and more general structures like the so-called S-and Z-sums. For our purpose some special cases of the algorithms A and B of this work are needed.
If an addend is not of the simple form S(a)/a m , but contains an offset in the denominator, which must not be symbolic, the summand can be reduced recursively This is a special case of Algorithm A. Of course, this is also valid for sums of depth greater than 2.
Algorithm B describes the recursive simplification of sums where products of harmonic sums appear.
can be rewritten by inserting the definition of S l . Then shifting variables, reordering of summations and partial fractioning yields to expressions which can be evaluated immediately and such of the same form as the original expression to which the procedure is applied again. Thus the sum is reduced recursively to simpler expressions and finally evaluated.

Products of Harmonic Sums
In Ref. [32] algebraic relations between harmonic sums have been studied. For harmonic sums with identical arguments the permutation relation holds: This allows us to express the product of two harmonic sums in the following way: For products with higher harmonic sums one has relations like More relations of this kind can be found in [33].

Harmonic Sums with Double Argument
In the following also sums with double argument S k (2m) will appear. These cannot be handled with the algorithms described up to now. In Ref. [34] transformations of such objects to simpler sums have been discussed. From this work we need the relation which relates the sum over S(2k) to two sums over S(k), one of which is alternating. The opposite operation, the transformation of a sum with simple argument to one with double argument, can be done by means of a special case of the refinement formula [34] S m (n) = 2 m−1 (S m (2n) + S −m (2n)) .
Generally, we can write a sum of the form in the following way: The sum over all even numbers is replaced by the sum over all integers minus the sum over the odd numbers. This is achieved by subtracting an alternating sum and the division by 2.

Multiple Infinite Sums
Finally, we need some formulae to convert nested sums to multiple infinite sums or vice versa. Two important relations are given by Eq. (48) can be obtained by splitting the sum on the left-hand side into two sums, which is visualized in Fig. 3, and paramaterizing the summation variable m = 2µ and m = 2µ + 1, respectively. This leads to sums of the type If the argument of the appearing harmonic sums is a sum of two or more variables, the following relations are helpful:

Convergence Acceleration
If the sums cannot be completely simplified analytically, one can at least get a numerical result by summing up the terms with the computer. However, in our case the sum converges very slowly. In this Section we briefly sketch the method of nonlinear sequence transformations, by means of which one can improve the convergence of a series. For further information we refer to Ref. [35] and references cited therein.
We consider a sequence {a n } and the corresponding infinite series with the partial sums which themselves form a sequence {S n } that converges to the value S. The question is now, whether a different sequence {S ′ n } (or {a ′ n }) exists that converges to the same value S but has a better convergence behaviour.
The sequence {S ′ n } is said to converge faster than {S n } to the common limit S, if In practice this means that S ′ N is a better approximation to the exact value of the limit than S N , if the partial sum of order N is known.
Convergence acceleration methods can be found with the help of model sequences. This means, one investigates a particular class of sequences and constructs a transformation, i.e. a prescription how to calculate the elements of a new sequence {S ′ n } out of the elements of {S n }, which is then exactly valid for the considered sequences. The prescriptions one has found in this way are often able to improve the convergence of different sequences of similar structure.
One of the oldest and best known methods is Aitken's ∆ 2 -process. The partial sums are transformed according to the rule ∆ is the forward difference operator, which is defined as [35,31,27] ∆f (n) = f (n + 1) − f (n).
Powers of ∆ are to be understood as multiple application of ∆, e.g.
Another method which can be found in the literature on numerical methods is Wynn's epsilon algorithm, which is defined by the recursion scheme We follow the notation of Ref. [35]. The lower index indicates the transformation order, the upper index labels the sequence elements. The difference operator acts only on the upper index.
Similar to the epsilon algorithm is Wynn's rho algorithm: where the interpolation points x k have to fulfill 0 < x i < x j for i < j and the sequence {x n } diverges to infinity. A possible choice is x n = n + 1.
An example, where the rho algorithm is very powerful, is the series defining When one sums the series up to the order i = 2000 one gets only 4 correct decimal digits. Application of the rho algorithm delivers more than 250 correct digits.
For the Levin transformations remainder estimations have to be taken into account. The Levin transformation ist exact for sequences of the form where S is the limit of the sequence and β > 0 is a parameter which can, in principle, be chosen arbitrarily and is usually set equal to 1. ω n are the remainder estimates. They are, of course, generally not known exactly, because the actual limit is unknown. Therefore one uses assumptions for the ω n which are customized to the given problem and which define different kinds of Levin transformations.
As an example, the ansatz ω n = a n a n+1 a n − a n+1 , n ∈ AE 0 , An enhancement of the epsilon algorithm is Brezinski's theta algorithm: It is now suggestive to apply such a procedure, that improves the convergence of a series, succesively. In Ref. [35] algorithms are discussed which allow the calculation of high transformation orders efficiently and with low memory costs.
The principle is to calculate for each new element of the series all possible elements of the matrix ϑ j i . So one always gets the highest possible transformation order. Details can be found in Ref. [35], where also many more acceleration methods are discussed.
For this work the algorithms mentioned above were programmed in FORTRAN. To achieve the desired numerical precision we used Bailey's library MPFUN [36,37], which allows calculations with arbitrary precision.

Numerical Result for the Nonplanar Topology
In this Section we sketch the application of the described methods to the nonplanar diagram N and show how a numerical result is obtained.

Coefficient of ε 1
The coefficient of ε 1 is a triple sum of rational expressions, where no harmonic sums occur. The sum can be evaluated analytically, which we sketch briefly in this subsection. The sum over l can be performed easily and we are left with a double sum. In the summands we get harmonic sums with arguments m and m + n.
The resulting rational terms and terms with S(m) can now be summed over n. The terms containing S(m + n) can be cast into an expression that can be summed over m by means of a transformation of variables. After the summation we are left with a simple sum with increased number of terms. The summands now consist of rational expressions, harmonic sums with simple argument S k (i + n), such with double argument S k (i + 2n) and products of the latter.
Again, the rational expressions can be summed up analytically by means of Mathematica or SUMMER. Expressions containing harmonic sums with simple arguments can be summed with SUMMER.
Expressions with harmonic sums with double argument can be transformed to terms calculable with SUMMER by means of formula (43).
By means of Eq. (44) the harmonic sum is converted to sums with double argument Then the sum is rewritten according to Eq. (46) and can now be calculated.
Finally, there are the products which can be evaluated with the help of Eq. (44).
The problem for the ε 1 coefficient can thus be solved completely and the result readsÑ which agrees with [10].

Coefficient of ε 2
In the expression for the ε 2 coefficient the terms in the triple sum contain the harmonic sums S 1 (n − l), S 1 (m), S 1 (l) and S 1 (m + n). All expressions can be summed over a variable which does not appear as an argument in S.

The result is then a double sum
where terms with harmonic sums of the following types occur: In the case of S i (m) and S i (m) S j (m) we sum now over the variable n which is not the argument of S. Expressions containing S i (m + n), S i (m + n) S j (m + n) can be handled by means of variable transformations as described in Section 3.7. The terms with S i (m) S j (n), S i (m) S j (m + n) can in the most cases be cast into a form to which algorithm B of Ref. [26] can be applied.
Thus the second summation yields a simple sum. Now the addends contain simple harmonic sums and products of two or three harmonic sums, some of which have double arguments: These terms can in principle be calculated by means of the described algorithms.
However, in some of the terms there are certain combinations of denominators and harmonic sums in which the summation variable is doubled. These can in some cases be be handled with the known algorithms, however, not in the general case. Therefore the problem of the order ε 2 cannot be solved analytically at present. To achieve this, more examinations in the subject of the summations are necessary.

Numerical evaluation
Since up to now an analytical solution is not possible with our method we try to find a numerical solution. To do this, we simplify the terms as far as possible and calculate the sums numerically.
For the numerical calculation the terms have been simplified with Mathematica and output in FORTRAN syntax. This code has been optimized by means of Maple and built in a FORTRAN program. To speed up the calculation this program was parallelized by means of MPI [38,39]. Finally, the program was translated with Bailey's program TRANSMP [37] to allow arbitrary numerical precision by means of the MPFUN library [36]. The programs store all partial sums.
To the sequence of partial sums we apply different convergence acceleration methods.
A different approach to a numerical solution, is the following: We begin with the representation of the sum (34) For a fixed n = ν the expression describes ν + 1 infinite sums Sums of that kind can easily be calculated with SUMMER. So the obvious procedure is to evaluate all contributions S n,l for all values of n up to an upper limit n max .
Thus we get a sequence of expressions consisting of rational numbers and zeta functions. These can be evaluated numerically to achieve a sequence of decimal To this we again apply the convergence acceleration algorithms. The first elements of this sequence are This procedure can be implemented easily in FORM or its parallelized version PARFORM [40,41].
We show in Table 1 Table 1 Numerical results for the ε 2 coefficient of N One can see, that the results of the Aitken and the Rho method agree in 21 decimal digits. These algorithms delivered the best result ‡ for the numerical calculations of the planar topology, where the exact result was known and could be used to check the results. Since the structure of the terms is similar, one can expect, that they are also suitable for the nonplanar diagram. The result of the Epsilon procedure agrees in 20 digits with the former. All results agree in 16 digits. In a very pessimistic interpretation at least these digits can be considered reliable.
The results of the numerical evaluation by means of the FORM-method agree with these numbers in 14 or 20 digits, depending on the acceleration method used.

Integer Relation Detection
Finally, we want to get an idea of the analytical result by means of "experimental mathematics". Looking at the known result for the planar diagram, we assume, that the result for N be also a linear combination of zeta values with integer coefficients, which have to be determined. This problem is known as integer relation detection [42].
In general a vector x = (x 1 , x 2 , ..., x n ) with real or complex elements is said to possess an integer relation, if integers a i exist such that at least one of them is different from zero and a 1 x 1 + a 2 x 2 + . . . + a n x n = 0 (75) ‡ Summing up 2000 terms we get 3 correct digits for P . Application of the rho algorithm yields 13 correct digits in this case.
An algorithm which can solve this problem is the PSLQ algorithm [43,44,45].
The input for this algorithm are in our case the numerically calculated number and numerical values for the constants, i.e. the zeta values. Then the program checks, whether there is an integer relation between the numbers within the desired accuracy. We used an implementation of this method by O. Veretin [46].
Assuming that the transcendentals are the same as for the planar diagram, as it is the case in the lower orders in ε, we find the following result for the 2nd order coefficient of the ε expansion of N: This agrees with the analytical result which was found in the meantime [19,20].
Comparison of this result with Table 1 shows, that all the digits printed for the rho algorithm are correct. This proves the power of this algorithm for our problem and justifies the confidence in the convergence acceleration techniques retrospectively.

Analytical Result for the Planar Topology
In this Section we show that the coefficient of ε 2 for the planar diagram can be found analytically by means of the described method. The steps we take are the same as described in Section 5. The structure of the sums and of the terms are similar and so is the procedure of the calculation. The difference is, the terms that cause problems in the case of N, because we do not know algorithms for their evaluation, do not appear in the case of P .

Conclusion
We described a method for the calculation of massless, dimensionally regularized Feynman integrals in their ε-expansion. It is based on the expansion of the integrand in terms of Gegenbauer polynomials. This yields sums over Gamma functions which can be expanded in harmonic sums. As a consequence one arrives at sums over rational functions and products of such with harmonic sums which can be simplified by the introduced summation algorithms and then be calculated numerically or analytically.
The procedure in principle is independent of the Gegenbauer polynomial technique. It can be useful wherever there are sums of the described form. The method should be expandable to more complicated problems, e.g. higher orders of the diagrams under consideration or more complex diagrams.
In the case of the nonplanar diagram we could not get to a complete analytical solution. To achieve this, further knowledge in the summation techniques is required.
For this work we combined the abilities of different computer programs, in particular the pattern matching facilities of Mathematica, the code generation functions of Maple, and FORM's power of managing very large expressions.
It is the idea to implement the method in one programming language to make it stand-alone. This could either be done in FORM or in C++. The summation algorithms of Ref. [26] have already been implemented in both, the FORM version has recently been published [50]. In [51] also a Mathematica implementation was introduced.
Finally the author wishes to thank M. Steinhauser and K. Chetyrkin for reading the manuscript and many useful comments.
This work was supported by the Graduiertenkolleg "Hochenergiephysik und Teilchenastrophysik" and by the Sonderforschungsbereich TR 9, "Computational Particle Physics".

A Properties of the Gegenbauer Polynomials
The Gegenbauer polynomials, also called ultraspherical polynomials, are the coefficients in the power series exansion of the generating function Information about these polynomials can be found in the books by Erdélyi et al. [47], Vilenkin [48] or Szegö [49]. The notation in this article follows Chetyrkin et al. [12].