The Method of Arbitrarily Large Moments to Calculate Single Scale Processes in Quantum Field Theory

We device a new method to calculate a large number of Mellin moments of single scale quantities using the systems of differential and/or difference equations obtained by integration-by-parts identities between the corresponding Feynman integrals of loop corrections to physical quantities. These scalar quantities have a much simpler mathematical structure than the complete quantity. A sufficiently large set of moments may even allow the analytic reconstruction of the whole quantity considered, holding in case of first order factorizing systems. In any case, one may derive highly precise numerical representations in general using this method, which is otherwise completely analytic.


INTRODUCTION
Single scale higher order QED and QCD calculations in the massless [1,2] and massive [3] cases at fixed Mellin moment n are given by polynomials of rational numbers and a series of a few special constants, the multiple ζvalues, and possible generalizations thereof [4]. This is irrespectively the case, whether or not the functional representation for general values of n obeys an equation factorizing in first order, [1], or not [5]. If one has access [6] to an algorithm, through which a large number of moments, e.g. N = 2000 or larger, can be calculated, the method of guessing, cf. [8], for holonomic problems, which often appear in physics applications, allows one to gain one large difference equation describing the corresponding problem [9]. If this equation is solvable in difference ring theory [10] one finds the solution for general values of n for this quantity without any further assumptions, e.g. made in [2,11]. In any case, significantly more moments will allow one to constrain the considered quantity much better numerically using approximation methods, e.g. through Chebyshev-polynomials or other interpolation methods.
In the following we describe an algorithm through which the system of differential equations, or associated to it, that of difference equations, available by the integration-by-parts relations [12], can be used to compute a large number of Mellin moments for the master integrals, and through them for the whole problem. The solution of the associated difference equations will need a relatively low number of initial values which have to be provided. The corresponding sequences of rational numbers mentioned above actually form the problematic part in gaining the general n result from the moments, since very involved, and in some cases yet unknown, functions span the corresponding sequences. In any case, the method allows either to find the one-dimensional distribution from a large but finite amount of moments, or at least to constrain it numerically at high accuracy.
We will give a brief illustration of the algorithm in case of a system of massive 3-loop master integrals, for which first a larger number of moments is generated, corresponding difference equations are found and solved for general values of the Mellin variable n.

THE ALGORITHM
Single scale master integrals can be represented as ana- We aim at computing a large number of coefficients I i (0), I i (1), . . . , I i (s). Usually the coefficients I i (n) depend on the dimensional parameter ε which itself can be expanded in a Laurent series in ε of a certain order o ∈ Z. We are interested in calculating the coefficients I (k) i (n) ∈ R of the expansionŝ up to a certain degree in ε, t i ∈ Z. More precisely, we want to compute for 1 ≤ i ≤ m the initial values for the ε-orders o ≤ k ≤ t i . In our approach we rely on the property that these unknown functionsÎ i (x) are usually described by a coupled system of first oder linear differential equations . . .
with D x = d/dx and where A is an m×m matrix with entries consisting of polynomials in ε and x and the entrieŝ r i can be given in form of the expansionŝ Another important assumption is that the coefficients r (k) i (n) can be determined efficiently by using either the method under consideration in a recursive fashion or by using, e.g., symbolic summation and integration techniques [10,13]. Further, we assume that for reasonable small numbers s ′ i the coefficients I i as a preprocessing step. Given this input we introduce the following efficient algorithm that computes the coefficients (2) up to large values of n.
(1) Using decoupling algorithms [14,15] transform the system (3) symbolically to one scalar linear differential equation where the a i (x, ε) and r i (x) are polynomials in x and ε; in addition, one obtains identities of the form for 2 ≤ i ≤ n where the e i,j,k (x, ε) and f i,j,k (x, ε) are rational functions in x and ε.
The algorithm proceeds now as follows: Compute in steps (2)-(4) the initial values forÎ 1 (x) by using the scalar differential equation (5), and compute afterwards in step (5) the initial values for the remaining integrals with the given formulas (6). (5) and eliminate D x by using the property D x ∞ n=0 h(n)x n = ∞ n=1 n h(n)x n−1 for a power series ∞ n=0 h(n)x n . Then by coefficient comparison w.r.t. x n and using appropriate shifts one gets a linear recurrence of the form d k=0 b k (n, ε)I 1 (n + k) = ρ(n, ε), with ρ(n, ε) = m j=1 l k=0 g j,k (n, ε)r j (n + k) for some l ∈ N where the b k (n, ε) and g j,k (n, ε) are polynomials in n and ε. Finally, divide the equation by a factor ε u for some u ≥ 0 in order to obtain updated polynomials b k (n, ε) in n, ε where not all b k (n, 0) with 0 ≤ k ≤ d are zero; the g j,k (n, ε) are now polynomials in n and Laurent polynomials in ε.
(3) We write the right hand side of (7) in the expanded representation with (4) can be computed efficiently by assumption, the coefficients ρ (k) (n) for sufficiently large k ∈ Z and n ∈ {0, 1, . . . , s} can be obtained explicitly without any cost.
(4) We proceed as follows; compare [16], Lemma 1. Plugging into (7) and doing coefficient comparison w.r.t. ε o yield the constraint Then we can compute with the first values I for n ′ = n − d ′ , n ≥ d ′ + δ all the other values in linear time. Now insert (9) with the explicitly computed values I 1 (n) with 0 ≤ n ≤ s into (7) and move these given values to the right hand side. This yields where the ρ ′(k) (n) for 0 ≤ n ≤ s and sufficiently large k are given explicitly. Now we are in the position to repeat this tactic iteratively to compute the remaining coefficients: next I o+1 1 (n) with 0 ≤ n ≤ s, afterwards I o+2 1 (n) with 0 ≤ n ≤ s, and eventually I ti 1 (n) with 0 ≤ n ≤ s.
(5) Now expand the rational functions f i,j,k (x, ε) and g i,j,k (x, ε) in (6) in a Laurent-series expansion w.r.t. ε and compute their coefficients as power series in x up to the necessary orders. Finally, combine all the explicitly given (finite) expansions in (6) using component-wise addition and Cauchy-products which yield the coefficients I The following remarks are in order; related considerations have been applied also to our algorithms to solve coupled systems in terms of nested sums over hypergeometric products [13,[17][18][19].
(i) Our algorithm requires that sufficiently many initial values/coefficients ofr i (x) andÎ 1 (x) up to the right ε-order are computed as a preprocessing step. The necessary bounds for these numbers can be determined by analyzing the formulas (7) and (6) accordingly.
(ii) For simplicity, we assumed that only one scalar differential equation arises, as it happens in most examples. In general, several scalar differential equations might arise, i.e., steps (2)-(4) have to be executed several times.
(iii) One can choose any λ with 1 ≤ λ ≤ m to determine a scalar differential equation inÎ λ (x). Different choices of λ might lead to different recurrences (7) with different orders d and formulae (6) of different size. In our calculations we analyze for each λ (1 ≤ λ ≤ m) the obtained symbolic formulae and choose that λ for the calculation steps (3)-(5) that serves us best, e.g., to minimize the required ε-orders for ther i (x) or minimize the the recurrence order d.
(iv) Often it is a challenge to determine the first initial values of I k i (n) to activate the recurrence formula (10). Thus it is highly desirable to keep the order d in (7) as small as possible. Since d gets smaller if the degrees of the a i (x, ε) w.r.t. x in (5) can be made smaller, the following refinement can be applied. Divide (5) by the polynomial h = gcd x (a 0 , . . . , a n ) in x. In most applications this reduces the degrees of the polynomials a i (x, ε) w.r.t. x heavily and thus produces recurrences with much lower order, e.g., d = 4 instead of d = 20 in (7). The price to be paid is to determine the coefficients ρ (k) (n) in (8) by extracting the nth coefficients from the right hand side of (5). Here the more involved calculation steps similarly as given in step (5) must be carried out.
All these aspects have been implemented within our package SolveCoupledSystem [13,18,19] which uses subroutines of the summation package Sigma [10] and the uncoupling package OreSys [15].

AN EXAMPLE
Let us consider a set of three 3-loop master integrals {Ĵ 1 (x, ε),Ĵ 2 (x, ε),Ĵ 3 (x, ε)} which obey a 3 × 3 inhomogeneous linear coupled system, cf. (3), as an illustrative example. They contribute to the massive operator matrix element A (3) gg [20], and are given in Eqs. (3.62-3.64). In a first step, we will use our above method to compute a large series of fixed moments for J 1,2,3 (n, ε), and in a second step we will use this information to guess recurrences for these functions. In our example the integrals J i have poles up to ε −3 and we would like to compute the moments up to the constant term O(ε 0 ). Since the number of moments n 0 needed to guess a recurrence is not known a priori, the moments are computed up to a reasonable number. If this is not sufficient, further moments can be computed reusing the previous calculations. E.g., using our algorithm the time to generate n 0 = 500 moments amounts to 48 sec and n 0 = 2000 moments require 569 sec in the present example. In a nutshell, higher moments can therefore easily be produced.
The guessing method to find the minimal difference equation requires the numbers n 0 , depending on ε and whether the respective set corresponds to the purely rational term or a respective contribution ∝ ζ i . The values are summarized in Table I. Here the highest number turns out to be n 0 = 89. The guessing method [8] yields a difference equation  Table II. The number of moments always includes a sufficiently large safety mar-gin to validate the corresponding difference equations, which may even be enlarged. We generate additional moments according to their order to provide the needed initial values to solve the difference equations.
In the present case, using the algorithm of Ref. [13], all difference equations can be solved completely using difference ring theory [10]. Besides rational terms in n the result for J 1 (n) and J 2 (n) is spanned by the harmonic sums [21] up to S 1,2 (n). J 3 (n), furthermore, contains the finite binomial sum multiplied by a rational term in n, see also [22]. All appearing letters building the nested sums forming the master integrals are found automatically. The above example is a simpler one. Let us mention that, for comparison, the computational time for the complete pole terms of more involved massive 3-loop examples requires the knowledge of about ∼ 1000 even moments [9] and is estimated to amount to several weeks [23]. In case the difference equations would turn out not to be first order factorizable using the algorithm of Ref. [13], at least the first order factorizable terms would be gained analytically, leaving the non-factorizable part behind for further analysis using different methods, like those for elliptic integrals and their potential generalization.
Having the above number of moments, Eq. (1) yields a first numerical approximation ofĴ k (x) also up to the required power in ε.
The present method is suitable to obtain precise numerical representations in case of various current massless and massive calculations at 3-and 4-loop order. In case one may solve the associated difference equations algebraically, one will even obtain the complete analytic result without any prejudice, i.e. a sufficiently large set of scalar moments allows one to fully reconstruct a one-dimensional distribution. This automatic method is therefore suited to carry out many more higher loop calculations in contemporary elementary particle physics in a completely or at least widely analytic form.