Calculating Massive 3-loop Graphs for Operator Matrix Elements by the Method of Hyperlogarithms

We calculate convergent 3-loop Feynman diagrams containing a single massive loop equipped with twist $\tau =2$ local operator insertions corresponding to spin $N$. They contribute to the massive operator matrix elements in QCD describing the massive Wilson coefficients for deep-inelastic scattering at large virtualities. Diagrams of this kind can be computed using an extended version to the method of hyperlogarithms, originally being designed for massless Feynman diagrams without operators. The method is applied to Benz- and $V$-type graphs, belonging to the genuine 3-loop topologies. In case of the $V$-type graphs with five massive propagators new types of nested sums and iterated integrals emerge. The sums are given in terms of finite binomially and inverse binomially weighted generalized cyclotomic sums, while the 1-dimensionally iterated integrals are based on a set of $\sim 30$ square-root valued letters. We also derive the asymptotic representations of the nested sums and present the solution for $N \in \mathbb{C}$. Integrals with a power-like divergence in $N$--space $\propto a^N, a \in \mathbb{R}, a>1,$ for large values of $N$ emerge. They still possess a representation in $x$--space, which is given in terms of root-valued iterated integrals in the present case. The method of hyperlogarithms is also used to calculate higher moments for crossed box graphs with different operator insertions.


Introduction
Massive on-shell operator matrix elements (OMEs) occur in the calculation of the Wilson coefficients in deeply-inelastic scattering, describing these quantities at large enough virtualities Q 2 m 2 together with the massless Wilson coefficients [1]. These OMEs are loop corrections to local composite operators being placed in graphs with massless external lines, which are on-shell. Their scale is set by the mass of an internal closed fermion line. Starting at 3-loop order, graphs with more than a single mass contribute [2,3]. The scale Q 2 from which on the asymptotic representation was found to apply at 2-loop order at the 1% level for the structure function F 2 (x, Q 2 ) is Q 2 /m 2 > ∼ 10, with m the heavy quark mass, cf. [1]. Here the asymptotic result was compared to the complete one [4] also containing non-universal power corrections. For F 2 (x, Q 2 ) this is a very acceptable kinematic range at HERA in case of m = m charm since at lower virtualities Q 2 < ∼ 20 GeV 2 still significant higher twist terms contribute [5][6][7].
Beyond NLO all massive OMEs have been calculated for a series of moments N = 10, (12,14) in the single mass case [8,9] for F 2 , F L and transversity, and the moments N = 2, 4, 6 for the contributions with two different masses [2,3,10] for F 2 at NNLO. With these results also all contributions to the unpolarized 3-loop anomalous dimensions ∝ T F were calculated independently for these moments and confirmed earlier results, cf. [11].
In case of the massive OMEs and Wilson coefficients at general values of the Mellin variable N all logarithmic contributions are available [12], to which also the 2-loop terms [1,13] up to O(ε) [14] contribute. 2 All O(T 2 F N F ) contributions were computed in [16,17]. This includes the two complete massive 3-loop OMEs A were calculated [18]. There are first results on the T 2 F -terms in the equal mass case [2,19,20]. In the polarized case the massive OMEs were computed to 2-loop order in Refs. [21,22]. In the calculation of these diagram classes the Feynman parameter integrals are reduced to multiply nested finite and infinite sums [23,24], using representations through hypergeometric functions and their generalizations [25] and Mellin-Barnes representations [26]. The sums obtained are then calculated using the packages Sigma [27], EvaluateMultiSums and SumProduction [28], applying also the algebraic and structural properties of harmonic sums [23,[29][30][31][32], their associated polylogarithms [33] , and special constants [34], including extensions to the cyclotomic [35] and generalized harmonic sum case [36,37]. These relations are encoded in the package HarmonicSums, cf. [37][38][39]. 3 Beyond the above topologies at the 3-loop level also ladder and Benz-type 4 , V -type and crossed box graphs contribute. In Ref. [41] we calculated diagrams of the 3-loop ladder topology of up to six massive propagators, including the most demanding cases. Not all of these graphs could be calculated using the above technologies.
In case the corresponding graph exhibits no poles in the dimensional parameter ε = D − 4, the method of hyperlogarithms has been devised for massless 2-point topologies with an off-shell external momentum in scalar field theory in Ref. [42]. This method allows to transform the Feynman-parameter integrals into special numbers, which are for the first loop orders linear combinations of multiple zeta values [34] and are given in terms of hyperlogarithms at unity argument. The integration is organized as a consecutive mapping into hyperlogarithms due to the linear structure of Feynman-parameters in these integrals within every integration step, which is being kept during the integration process. In the present paper we generalize this method allowing for local operator insertions. Furthermore, we consider the case of massive diagrams in which a higher nesting of Feynman-parameters is generally expected if compared to the massless case. I.e. the formalism may lead to structures beyond linearity at an earlier stage than in the massless case. The local operator insertions introduce a new degree of freedom, the Mellin variable N . The corresponding Feynman diagrams are given in terms of sum-representations. In most simple cases harmonic sums emerge. More involved cases lead to generalized sums over rational alphabets, and also nested cyclotomic and binomial sums, as will be shown below. Interestingly, for fixed integer values of N the corresponding graphs evaluate to rational numbers, weighted by multiple zeta values for the loop-level considered in this paper, similar to the case in the original approach [42]. One may calculate moments up to N = 9 even for the most complicated 3-loop graphs which emerge in the present physics project. These moments can be checked by very different methods based on the codes MATAD [43] and qexp [44] at lower values of N . The method works, since the numerator functions are polynomials in the Feynman parameters at fixed values of N . Partial fractioning may be performed until one obtains denominator functions only. However, the above number of moments is usually still far too low to try the reconstruction of the general N behaviour using the method described in [45].
Introducing an auxiliary parameter x, the local operator insertions may, however, be resummed such that a generating function is obtained, which is expressed in terms of hyperlogarithms L a (x). In turn the N th Taylor-coefficient of this function has to be obtained analytically. The last step can be performed in some cases using HarmonicSums directly. In more complex situations associated difference equations of larger order have to be established and solved using Sigma [27].
The paper is organized as follows. We describe the extension of the method [42] to massive operator matrix elements at 3-loop order in the presence of local operator insertions in Section 2. The method is applied to convergent Benz-type and related diagrams in Section 3, also discussing practical aspects. Here we also derive the asymptotic representations of the individual graphs, which is necessary for their representation for complex values of N needed to perform the Mellin inversion in practical applications [23,46]. In Section 4 we calculate graphs of the 3-loop V -topology with five massive propagators. They may be considered to emerge from either a ladder-or the crossed box-topology by removing one line. While in the former case conventional structures are obtained, in the latter case new nested sum-types emerge, which contain weights due to binomials of the type 2i i both in the numerator and denominator. In the calculation root-valued structures in the auxiliary parameter occur in the last step which are responsible for these new hypergeometric terms. Aspects of the Mellin-inversion of the contributions from binomially weighted nested sums are discussed in Section 5. In Section 6 we apply the algorithm to calculate three crossed box-topologies for fixed integer values of N to demonstrate the applicability of the present algorithm also for these diagrams. Section 7 contains the conclusions.

The Formalism
We consider massive Feynman diagrams at l = 3 loops with operator insertions in D = 4 + ε dimensions. One may represent the Feynman parameter integral I G of a graph G in terms of Schwinger parameters using Symanzik [47] or Kirchhoff polynomials [48], cf. [49]. The Feynman rules are given in [8,50], including those for the operator insertions. The integral is given by : Here a i denote the powers of the different propagators, a = i∈edges a i . According to the Cheng-Wu theorem [51] the sum of Schwinger parameters over an arbitrary subset of edges E in G may be set equal to one, as expressed by the δ-distribution in (2.1). We associate to the graph G the graphG which is obtained by closing the external lines. While M G is given by the sum of all Schwinger parameters which are attached to a massive line, the graph polynomial Ψ G and the operator insertion OP i (α i , N ) obey the following graph theoretical descriptions. For a graph with n v vertices and n e edges we define the n e × n v graph incidence matrix if the edge e is not connected to vertex v .
We choose ε G as the matrix n e × (n v − 1)-matrix obtained from (2.2) by removing one arbitrary column. ε G is thus not uniquely defined and depends on the direction of the edges and the choice of the removed column. The graph matrix M G reads 3) The first graph polynomial Ψ G is given by Ψ G = − det(M G ). Although the matrix M G is not uniquely defined Ψ G , is independent of the possible choices for M G . If I,J,K are sets of edges in the graph G and I and J are of equal length, |I| = |J|, we define the Dodgson polynomials by with M G (I, J) being the matrix M G after removing all rows corresponding to the edges in I and all columns corresponding to the edges in J. If K is empty we omit it and write Ψ I,J G . The different operator insertions used in the present paper are expressed in terms of Dodgson polynomials given in Figure 1 for the examples studied in the present paper. The Dodgson polynomials Ψ I,J G,K are only defined up to a sign, which generally depends on the orientation of the edges in ε G and also on the column which has been removed to define M G . For the present paper we were able to choose Ψ I,J G,K = det M G (I, J) αe=0 ∀ e∈K if the directions of the edges correspond to the Feynman rules of Refs. [8,50].
Under certain conditions, Feynman parameter integrals, being convergent in D = 4 dimensions, can be cast into a linear combination of hyperlogaritms L( a, z) [55][56][57]. In the following we will outline the corresponding formalism, extending the algorithm [42], given originally for massless Feynman diagrams to those with also massive lines and local operator insertions.
Let σ be a set of distinct points in C and A = {a 0 , a 1 , ..., a N } an alphabet. We form words described by a out of the elements of A, where each letter corresponds to an element in σ. The elements in σ may be constants or rational functions of further parameters. The hyperlogarithms are defined by Here {...} denotes an ordered set. The weight w of a hyperlogarithm is given by the number of letters in a. The hyperlogarithms satisfy shuffle relations, cf. e.g. [31], In the shuffled index set one sums over all hyperlogarithms with indices such that the relative order of the indices in a 1 and a 2 is preserved. An example is given by The derivatives w.r.t. the argument z is The general tactic is to treat the inner most integral first and to transform the integrals from inside to outside in terms of hyperlogarithms. Let us start with the inner most integral and turn to the construction of the antiderivative (primitive functions) of products of rational functions R(z) = N (z)/D(z) and hyperlogarithms. Here we will assume that D(z) factors linearly, i.e. D(z) = k (z − a k ) l k , l k ∈ N. If there exists one integration order for a graph G for which this property is found in each integration step such a graph is called to be linear reducible. The consecutive decomposition of the multiple integral into a sequence of these steps is called Fubini sequence. Whether or not this decomposition exists can be checked a priori with reduction algorithms given in Refs. [42,53] by which also the requested order of integration is delivered. Applying the shuffle relation and partial fractioning one arrives at expressions of the form For n = −1 again the hyperlogarithm L ({−b, a 1 , a}, x) is obtained. Otherwise one applies integration by parts where in the last term the weight of the hyperlogarithm is reduced by one. Applying this technique recursively, all integrals can be written in terms of hyperlogarithms that have to be evaluated at its integration bounds in the α-representation (i.e. at 0 and ∞). The challenge is now to perform this evaluation, more precisely to calculate the limits. To accomplish this task, we actually calculate the series expansion at 0 and at ∞ and express the result again in terms of hyperlogarithms afterwards. This finally enables one to apply the presented method for the next integral. Next we consider series expansions of the hyperlogarithms around z = 0 and for z → ∞. A hyperlogarithm of weight w satisfies series representations of the form Following [42] it is suitable to define the restricted regularization RReg z→{0,∞} given by the constant part of the generalized series expansion RReg z→0 L({a 1 , · · · , a n }, z) = c RReg z→∞ L({a 1 , · · · , a n }, z) = c The series expansions are constructed as follows. One first differentiates w.r.t. the argument of the hyperlogarithm and then performs the series expansion of the derivative, which is of lower weight. After this the antiderivative is calculated and the respective integration constants are fixed. We denote the series operator by Ser (k) y→∞ , up to terms of O y −k log w (y) . One obtains For example, one finds The same method is applied to construct the series representations for hyperlogarithms of higher weight. 5 We now line out how the integration constants can be transformed, which is necessary in the applications. Derivatives for a the variable t of which the letters a i (t) in the index-set of the hyperlogarithms may depend, are computed as follows : .

(2.23)
Note that taking the derivative with respect to the argument or an inner variable of the hyperlogarithm always yields expressions which contain only hyperlogarithms of a lower weight. To prepare the next integration step, the constants c (∞) 0,0 ({a 1 , · · · , a n }) = RReg y→∞ L ({a 1 , · · · , a n }, y) (2.24) have to be rewritten in terms of hyperlogarithms, such that the next integration variable does not appear in the respective index set. This is done by differentiating, rewriting the now weightreduced expression and then forming the antiderivative again. Let us consider the example Algorithms to obtain closed forms for these expansions are known and have been implemented into the computer algebra package HarmonicSums [37][38][39].
with ζ k = ∞ l=1 1/l k , k ∈ N, k ≥ 2 the Riemann ζ-function. Special care has to be taken when evaluating constants like c (∞) 0,0 (a 1 , · · · , a n ) which contain letters of the form x −i f (x) with f (x) = 0 as x → 0 or trailing letters of the form x i f (x), with lim x→0 f (x) being finite. In all other cases RReg x→0 L (a 1 , · · · , a n , y) is just obtained by taking the limit x → 0 under the integral. In the first case the limit x → 0 does not commute with y → ∞. If a hyperlogarithm does not have any trailing zero in its index set, we may substitute the integration variables z i → az i in (2.9) to obtain L ({a 1 , · · · , a n }, z) = L ({aa 1 , · · · , aa n }, az) . (2.28) In other cases trailing zeros have to be removed by means of the shuffle algebra first, e.g., By definition Ser (0) z→∞ L ({f 1 (x), · · · , x i f n (x)}, z) does depend on the variable z = yx i only logarithmically and the operation RReg y→∞ in the second last step is easily performed. In the case of trailing letters of the type x i f (x) with f (x) finite as x → 0, the limit x → 0 does not commute with the implicit limits contained in the definition of the hyperlogarithm. Here we apply the identity on the parts containing the respective letters as many times as needed. The identity (2.32) is derived by considering the change of integrations variables z i → z i x in (2.9). We illustrate this in the following example : The previous steps are repeated for all further integration variables until we have rewritten all constants in a way suitable for the following parametric integrations. That far we have described the algorithm for a finite loop diagram built of propagators and vertices for a renormalizable quantum field theory. The present application is more general as also local operator insertions shall be dealt with. A consistent set of Feynman rules in case of Quantum Chromodynamics has been presented in Ref. [8]. As a consequence of the light-cone expansion [58] the local operator insertions emerge as polynomials of degree N , N ∈ N, as has been outlined above. For any integer value the present formalism can be applied through which the moments of the corresponding OME are obtained. With growing values of N both the requested CPU time and memory to perform this computation will grow significantly, usually with a nearly constant factor by going from N → N + 2. All finite 3-loop topologies can be dealt with this method up to a certain moment, i.e. the present method is equivalent for finite diagrams to MATAD [43], which, however, can handle divergent graphs as well. In Section 6 we will illustrate this for the most complicated graphs in the present project.
To use the present method also in case of general values of the Mellin variable N , the following resummation into a generating function in the parameter t of the operator-polynomials is applied, cf. [41] : (2.34) Let us illustrate the derivation of the generating function for an operator insertion on a 3-vertex. It is of the structure The infinite resummation results into The generalization to the case of l-leg operator-insertions is straightforward. It leads to (l − 1)additional propagator terms, now containing also the variable t. In this way structures are obtained which are in a form suitable for the above algorithm. In case the auxiliary parameter t does not destroy linearity in the consecutive integration of Feynman parameters, finally a representation of the generating functions by hyperlogarithms L w (t) is obtained.
The following representations hold for the three different operators given in Figure 1 : The solution for the general Mellin variable N can finally be obtained by calculating the N th expansion coefficient of the generating function. This usually requires to solve associated difference equations. Respective algorithms are encoded in the packages Sigma [27], EvaluateMultiSums, SumProduction [28] and HarmonicSums [37][38][39]. We finally would like to note that for a fixed value of N all massive 3-loop QCD two-point topologies turned out to be linear reducible in the case of a single mass scale m. If we introduce generating functions this changes drastically. Some diagrams remain linear reducible, others can be transformed into linear reducible diagrams via a variable transformation. There are, however, also cases for which no sequence could be found to restore linear reducibility. One of the finite diagrams we would like to calculate is the scalar graph shown in Figure 2 in Section 4. For this diagram no completely linear reducible integration order exists a priori. 6 The linearization of some quadratic forms occurring can be performed introducing complex letters. A final quadratic form appears in the last step only and can be dealt with remapping the tracing variable to gain linear reducibility by .
The final expression will consist of hyperlogarithms in the new variable x = t/t + 4. More evolved techniques have to be applied to obtain the N -space representation, see Section 5.
We now turn to the calculation of specific finite 3-loop topologies applying the above methods.

Benz-Graphs
Let us first consider so-called Benz topologies. A first example is given in Figure 3. Here all powers of the propagators are chosen as ν i = 1. Using the method described in Section 2 one obtains the following expression : Here the global N -dependent factors stem from pre-manufacturing. The hyperlogarithms in (3.1) are even harmonic polylogarithms (HPLs) over the alphabet {0, 1, −1} [33]. Considering (3.1) as a power series in x, the N th coefficient of this expression in x has to be extracted analytically in order to recover the original integral. This can be achieved using the GetMoment function of the package HarmonicSums, cf. [37]. One may also use guessing-methods to obtain the corresponding difference equation based on a huge number of moments, cf. [45], and obtain the N th coefficient by solving this equation using Sigma [27]. The N th Taylor coefficient of (3.1) is given as the following representation in harmonic sums : The harmonic sums are denoted by [29,30] and we use the short-hand notation S a (N ) ≡ S a . For all finite sum structures one easily derives the recursive shift relation All harmonic sums can be written in terms of polynomial factors in S 1 (N ) and those [23], which have representations by factorial series [60]. The singularities of these sums are located at the non-positive integers, implying that these are meromorphic functions. Furthermore the physical expressions may exhibit singularities due to rational factors. The rightmost singularity is determined by the spin of the particles involved. In case of massless spin-1 (1/2, 0) particles singularities up to N = 1 (0, −1) can occur. The asymptotic representation of both types of sums can be uniquely determined and is automated by the code HarmonicSums. The asymptotic representation and the shift-relation (3.6) allow the analytic continuation of integrals like I 1 (N ) into the complex plane. The uniqueness of the analytic continuation can be proven by an extension of Carlson's theorem [37]. It is carried out either from the even or odd integers N in the sum expression, depending on the crossing relations of the process described, cf. [58]. Therefore alternating sums and factors (−1) N have a definite meaning prior to the analytic continuation N ∈ C. For Eq. (3.2) the asymptotic expansion is given by  Following the above algorithm, integral I 2 (N ) defined by the graph in Figure 4, yields : This integral contains generalized harmonic sums and also terms of the O(2 N ), which cancel in the asymptotic expansion.
(3.9) Diagram I 3 (N ) differs from diagram I 1 (N ) by moving the operator insertion to one propagator to the right. The result obtained is much more simple than for I 1 (N ), cf. (3.2), and is given in terms of a few harmonic sums only,   Integral I 4 (N ) is given by (3.14) Despite diagrams I 4 and I 5 are topologically quite similar, their result turns out to be structurally different. Integral I 5 (N ) is given by (3.16) Finally we consider diagram 6 as an example for convergent Benz-graphs. Applying the above algorithm one obtains : Here C 1 and C 2 are group-theoretic factors accounting e.g. for color degrees. We leave them unspecified since only scalar graphs are calculated, see Figure 1. The polynomials in (3.17) read (3.31) The asymptotic expansion of I 6 is given by I asy 6 (N ) = C 1 For all the above graphs, irrespectively of their concrete representation at integer values of N , which is of different complexity, the shift relation N → (N − 1) for N ∈ C can be established through simpler functions correspondingly, for which the analytic continuation has been worked out in Refs. [23,35,37]. In case of harmonic sums and cyclotomic harmonic sums the singularities are located at N ∈ Z, N < 1. The rational pre-factors may induce also singularities at N = 1. The generalized harmonic sums in I 2 , I 4 and I 6 have already been studied in (3.36-3.40) in [41] giving the corresponding Mellin representations. They partly appear together with the pre-factor 2 N . As has been seen above, the corresponding asymptotic representations of I 2 , I 4 and I 6 are well behaved. We still have to determine the positions of the poles of these sums in the complex plane. The following integrals have to be considered: The last integral in (3.33) is analytic in C for any finite range. Thus the singularities of S 1 (2; N ) are those of S 1 (N ); the exponential growth of the sum for N → ∞ is canceled by other terms in the integrals I 2,4,6 . The second integral in the sum has a factorial series representation [60]. Here Li n (x) = ∞ k=1 (x k /k n ), n ≥ 0 denotes the polylogaritm. The singularities are thus located at the non-positive integers. This also applies for the sum S 1,1,2 (2, 1/2, 1; N ), related to the integrals Here H a (x) denote the harmonic polylogarithms over the alphabet {0, 1, −1} [33]. Next we consider with S 1,1 (m) = [S 2 1 (m) + S 2 (m)]/2. The representations of the harmonic sums [30] imply that (3.36) converges absolutely, with poles at −(N + l) ∈ N\{0}.

V-type Diagrams with Five Massive Propagators
Another genuine 3-loop topology is represented by the V -type diagram shown in Figure 2. According to the Feynman rules given in Figure 1 it consists out of two contributions, which are labeled by the constants C 1 and C 2 . One may consider these terms as being obtained by (a) either expanding one line of a ladder graph or (b) the crossed box graph, cf. Figure 6c, by applying the light-cone expansion. In the α-representation these graphs are given by where and the different graph polynomials read (4.4) The integral I 7a , stemming from a former ladder-like topology, is expected to have a representation and complexity of other ladder-type diagrams considered in Ref. [41] before. We first obtain the representation in terms of hyperlogarithms : The generating function representation is given by harmonic polylogarithms only. From (4.5) the N th Taylor coefficient is derived using the GetMoment option of HarmonicSums. I 7a (N ) is represented in terms of harmonic sums up weight w = 5 : (4.6) The asymptotic representation of integral I 7a is given by and shows a regular behaviour. Integral I 7b (N ), related to crossed-box topologies by one additional propagator expansion, conversely leads to new structures.
First we derive the representation ofÎ 7b (x) (4.2) in terms of iterated integrals containing the auxiliary parameter x. We define The result is given by 1405 different hyperlogarithms. The corresponding expression is too long to be given in full form here. Instead we show a series of typical terms to illustrate different contributing functions : ... The index sets of the hyperlogarithms contain the letters 10) x or r as argument and reach weight w = 5.
In the last step of integration in determining (4.9) root-valued letters appear. Both due to the massive case studied here and the presence of the local operator insertion in the present case no complete Fubini sequence is obtained in the first place. However, transformation (4.8) establishes linear reducibility once agian and the corresponding integral can be solved.
To derive the N th Taylor coefficient from (4.9) has not been straightforward. Here we have chosen two ways. In a more simple approach we generated fixed Mellin moments from (4.9) and used the method of guessing [61] to derive a corresponding difference equation, cf. also [45]. We were able to generate 1500 moments. About 800 moments were finally needed to establish the difference equation. With Sigma [27] this difference equation could be solved in a time of 2000 seconds, through which the N th Taylor coefficient has been obtained. The method of guessing mostly delivers correct results with a failure estimated to be ∼ 10 −60 [61], yet it is not exact. Therefore we also derived from (4.9) the N th coefficient using Sigma [27] and HarmonicSums [37][38][39]. This computation requested two days of computation time confirming the result obtained by the method of guessing : The integral I 7b (N ), beyond the harmonic [29,30] and generalized harmonic sums [36,37] also contains a series of finite binomially and inverse-binomially nested sums, summing over generalized harmonic sums. These structures emerge from the hyperlogarithms containing the set of letters in the alphabet (4.10) beyond those of harmonic polylogaritms and the root-function r(x) in the argument. It is the strength of packages like Sigma [27] based on general summation algorithms operating on difference fields to find the new sum-structures. Furthermore, the representation (4.11) is given by sums being transcendental to each other. Here we made use of sum representations having been introduced previously in Refs. [29,30,[35][36][37] which occur at lower levels of the sum hierarchy implied by Feynman integrals.

Analytic Continuation of Binomially Weighted Nested Sums
To obtain the analytic continuation of the binomial sums as given in (4.11) we first derive their representation in terms of a Mellin transformation Individual nested sums then usually are given as a linear combination where the constants c j and functions f j (x) do not depend on N . The functions f j (x) are defined in terms of iterated integrals. As starting point we only need the following basic integral representations : where the letters f w i (x) are given by Here and in the following we refer to the notation of Ref. [62]. From the Mellin transforms (5.3-5.5) we can obtain integral representations for the nested sums step by step. In general the computation proceeds as follows. Starting from the innermost sum we move outwards maintaining an integral representation of the sub-expressions visited so far. For each intermediate sum this first involves setting up an integral representation for the summand a j (N ) of the form (5.2). This may require computation of Mellin convolutions, which we will describe in more detail below. Next we obtain an integral representation of the same form of by Mellin convolution with the result for the inner sums computed so far. Then by the summation property we obtain an integral representation for the sum (5.8). These steps are repeated until the outermost sum has been processed. Now, we take a closer look at how we compute Mellin convolutions, which is the most challenging part of the calculation. Formally, we rely on the convolution formulae which give us a definite integral depending on a continuous parameter and which can be written in the form In order to obtain a closed form for this integral, we first set up a differential equation satisfied by F (x) and then obtain a solution of this equation satisfying appropriate initial conditions. In the first step we exploit the principle of differentiation under the integral. If we have a relation for the integrand f (x, y) of the form for some coefficients c i (x) independent of y and some function g(x, y), then by applying the operator 1 x dy this gives rise to a linear ordinary differential equation for the integral F (x) x) + additional boundary terms. (5.14) Proper care has to be taken for evaluating the right hand side of this relation in the presence of singularities. There are several computer algebra algorithms for different types of integrands f (x, y) which, given f (x, y), compute relations of the form (5.13). They either utilize differential fields [63][64][65] or holonomic systems and Ore algebras [66][67][68]. For obtaining solutions to the generated differential equations the following two observations are crucial. All differential equations obtained during our computations factor completely into first-order equations with rational function coefficients and, moreover, these factors all have algebraic functions of degree at most two as their solutions. These two observations imply that solutions are of the form where r i (x) are rational functions and p i (x) are square-free polynomials. We define the iterated integrals H * w (x) by and f j (x) are the corresponding basic functions, which partly contain root-valued denominators. Using a dedicated rewrite procedure [69] based on integration by parts we can write a basis of the solution space in terms of the iterated integrals which is then used to match initial conditions. For the representation of integral I 7b 32 different letters f j (x) are needed, cf. [62]. As an example we consider the representation for the following double sum : Here the last Mellin-transform at argument N = 0 takes the value 2+(8/ √ 5)[ln( √ 5 − 1) − ln (2)], while the former one is a new constant, beyond the (cyclotomic) multiple zeta values. The letter f w 8 is given by To perform the asymptotic expansion of the Mellin-transforms we use the representation One may expand f (e −z ) e −z at z = 0 and integrate (5.18) term-wise to obtain the asymptotic expansion for |N | → ∞, arg(N ) = π using for c > −1 and k ∈ N. These expansions are automated in the package HarmonicSums. With these prerequisites at hand the asymptotic expansion of (4.11) can now be performed. It turns out, that part of the individual sums contributing to (4.11) diverge ∝ 8 N , 4 N and 2 N for large values of N . In case of the present scalar integral I 7b (N ) the terms ∝ 8 N and ∝ 4 N cancel, while some of the terms ∝ 2 N remain. We also have checked the principal divergence of this graph for N → ∞ numerically. In the physical case, accounting for all color and numerator structures, also these terms are expected to cancel between the different diagrams. Due to the contributing large class of new sums one expects also that a series of new constants beyond the multiple zeta values [34], generalized (cyclotomic) zeta values [35,37] contribute, see also [62].
The asymptotic representation of I 7b (N ) reads :   Here the constants A 1 and A 2 are given by We have expressed part of these constants numerically up to five generalized harmonic polylogarithms at x = 1. We checked using PSLQ [70] that no integer relation between these HPLs based on 100 digits is found. The numerical values of these constants can be derived from the following one-dimensional integral representations referring to classical polylogarithms x + 1 = −0.04281095672416394220 The numerical parts recruit from 20 one-and 17 two-dimensional integrals, which will be given in [62] in explicit form. One example reads Furthermore, beyond the usual multiple zeta values [34], also the following constants contribute Some of the latter constants express infinite cyclotomic harmonic sums [35] or represent the infinite sums

Moments for Crossed-Box Graphs
Using the method of hyperlogarithms also fixed moments of convergent graphs can be evaluated. The method relies on partial fractioning of the operator polynomial induced by the operator. Correspondingly, for large values of N , the number of terms grows exponentially. The calculation time and the requested storage are growing significantly. To illustrate the potential of the method in this respect we select the possibly most complicated graphs of the present physics project belonging to crossed box topologies.
While for more simple topologies more moments can be calculated given typical CPU times of various hours to days, in case of the above topologies the 9th moments could be calculated in about 8 hours requesting a storage of 35 Gbyte. The 10th moment would have needed storage of more than 200 Gbyte RAM due to the intense use of partial fractioning. Since the algorithm is implemented as Maple-code the available RAM is the limiting parameter, unlike the case e.g. for FORM-programs, using also fast external discs [72]. In comparison, the FORM-based program MATAD [43] allows to calculate a few higher moments as well having the same time-and storage resources.  Table 1: Moments of the finite crossed-box graphs (a-c) shown in Figure 9.

Conclusions
It has long been noticed that many results for zero-and single-scale processes in renormalizable Quantum Field Theories can be expressed in terms of iterated integrals or nested harmonic sums at the lower loop level [73]. Ideally, a direct method was sought for to arrive at these results right form the Feynman parameterization of the contributing diagrams. In case of convergent Feynman integrals the method of hyperlogarithms provides this way in case a Fubini sequence can be found for the diagram being considered. In the present paper we have extended this method to the case of massive diagrams including local operator insertions. The calculation of fixed moments does not pose a theoretical problem, since the expressions can be reduced in principle by applying partial fractioning. With growing values of N the complexity of the expressions rises significantly such that the corresponding number of terms cannot be swallowed even by modern computers anymore. To extend the present abilities, special software implementations outside coding systems based on Mathematica and/or Maple are necessary, to free the main storage and allow the use of fast discs to store intermediary results being processed subsequently.
For general values of the Mellin variable N at three-loop order in Quantum Chromodynamics topologies contribute also, for which root-valued letters appear in the alphabet. If these can be traded for the argument of the hyperlogarithm, the method remains applicable. This is, however, not the case for all massive 3-loop topologies. On the other hand, a remarkably wide class of diagrams can be calculated using the method of hyperlogarithms. At the technical side the operator insertions are mapped to propagator-type factors referring to a representation in terms of generating functions. At the end of the calculation the N th expansion coefficient has to be determined analytically for which techniques are available in the Mathematica-package HarmonicSums. In some of the graphs multiply nested sums weighted by binomials of the type 2i i in the numerator and denominator occur. To construct the analytic continuation of these sums to N ∈ C their asymptotic expansion for |N | → ∞, arg(N ) = π has to be calculated analytically. This requires the analytic Mellin-inversion of the corresponding sum expressions. We used Risch-algorithm methods to compute the corresponding iterated integrals, which request a larger amount of root-valued letters, cf. Ref. [62] for details. Also a series of new special constants beyond those of the multiple zeta values and their cyclotomic and generalized sum generalizations emerges in these expressions. Operating in difference fields and using the Rischalgorithm we arrive at minimal representations algebraically keeping only functions with relative transcendence to each other. The present methods also allow the representation of the integrals calculated in the present paper in x-space. Detailed transformation algorithms and results are given in [62].
The present analysis deals with convergent Feynman integrals only, while most of the Feynman graphs exhibit poles in the dimensional parameter ε = D − 4. The calculation also of these diagrams requires a suitable regularization to be carried out first and still needs a thorough algebraic implementation. The major limiting factor for a general application of the algorithm to massive problems, including local operator insertions, at present consists in the emergence of root-valued letters already at intermediate steps of the algorithm. A thorough mathematical treatment of these structures may be the subject of future investigations.