The O(\alpha_s^3 T_F^2) Contributions to the Gluonic Operator Matrix Element

The $O(\alpha_s^3 T_F^2 C_F (C_A))$ contributions to the transition matrix element $A_{gg,Q}$ relevant for the variable flavor number scheme at 3--loop order are calculated. The corresponding graphs contain two massive fermion lines of equal mass leading to terms given by inverse binomially weighted sums beyond the usual harmonic sums. In $x$-space two root-valued letters contribute in the iterated integrals in addition to those forming the harmonic polylogarithms. We outline technical details needed in the calculation of graphs of this type, which are as well of importance in the case of two different internal massive lines.


Introduction
The precision determinations of the strong coupling constant α s (M 2 Z ) [1] and the parton densities, cf. e.g. Refs. [2], in deep-inelastic scattering require the knowledge of the heavy flavor corrections to 3-loop order. The heavy flavor corrections were calculated at NLO in semi-analytic form in [3] 1 . To avoid contributions of higher twist, the analysis has to be restricted to large enough values of Q 2 . It has been shown in [5] that for Q 2 > ∼ 10 m 2 , with m the heavy quark mass, the heavy flavor contributions to the structure function F 2 (x, Q 2 ) are very precisely described using the asymptotic representation in which all power corrections ∝ (m 2 /Q 2 ) k , k ∈ N + are neglected. In this limit the heavy flavor Wilson coefficients can be calculated analytically. They are given by convolutions of massive operator matrix elements (OMEs) and the massless Wilson coefficients, cf. Ref. [5,6]. The massless Wilson coefficients are known to 3-loop order [7]. In the past the asymptotic O(α 2 s ) corrections were calculated in Refs. [5,[8][9][10][11][12][13] in the unpolarized and polarized case, including the O(α 2 s ε) contributions, and in [14] for transversity. The heavy flavor corrections for charged current reactions are available at 1-loop and in the asymptotic case at 2-loops [15].
At 3-loop order, a series of moments has been calculated for all massive OMEs for N = 2... 10(14) contributing in the fixed and variable flavor scheme [6]. All logarithmic terms to 3loop order including the contributions to the constant term due to renormalization have been computed in Ref. [16]. The 3-loop heavy flavor corrections to F L (x, Q 2 ) in the asymptotic case were calculated in [16,17]. First results for general values of N have been obtained for all OMEs for the color factor N F T 2 F C F,A [18,19] and 3-loop ladder, Benz-, and V -topologies [20,21]. First α 3 s T 2 F C F,A -contributions at general N were calculated for the flavor non-singlet and pure-singlet terms in [22] for two heavy quark lines carrying the same mass. Furthermore, the moments N = 2, 4, 6 in case of the OMEs contributing to the structure function F 2 (x, Q 2 ) with two different heavy quark masses were computed in [22,23]. In all the above cases the massive OMEs are calculated for external massless partons which are on-shell. Recently, the complete 3-loop OMEs A gq , A NS qq,Q and A PS Qq and the associated Wilson coefficients in the asymptotic region have been calculated in Refs. [24,25]. Also the case of massive on-shell external lines has been treated in [26] recently.
In the present paper we calculate the O(α 3 s T 2 F C F,A ) corrections to the massive OME A gg,Q with local operator insertions on the gluonic lines at general values of N. This matrix element is of importance to establish the variable flavor number scheme (VFNS) at 3-loop order. The terms of O(α 3 s T 2 F C F,A ) derive from graphs with two internal massive fermion lines of equal mass. Unlike the foregoing 3-loop results for massive OMEs at general values of N [16-18, 24, 25] new functions beyond the harmonic sums [27] appear, which belong to the finite nested binomially weighted harmonic sums [28]. Here they are of the type 2 (1. 1) which have been considered in [30] before. Here S a (N) denotes the nested harmonic sum More involved sums of this type contribute to the massive V -topologies, cf. Ref. [21]. For a larger class of diagrams the calculation of the corresponding graphs is performed using Mellin-Barnes representations and requires cyclotomic harmonic sums and polylogarithms in intermediary steps. The corresponding nested sums are then solved using the summation and representation techniques encoded in the packages Sigma [31], HarmonicSums [32][33][34], EvaluateMultiSums, SumProduction [35], and RhoSum [36]. For a few Feynman diagrams, it proved to be efficient to calculate them using integration-by-parts [37]. The corresponding master integrals were computed applying systems of linear differential equations. The paper is organized as follows. In Section 2 we discuss the structure of the gluonic operator matrix element. At O(α 3 s T 2 F C F,A ) 39 Feynman diagrams contribute. The calculation methods to obtain the result at general values of the Mellin variable N are outlined in Section 3 in detail. In Section 4 we present the results for the OME and also obtain the contributions ∝ T 2 F C F,A to the gluonic 3-loop anomalous dimension γ gg . Section 5 contains the conclusions. In the Appendix we present the results for a series of scalar integrals which emerge in the present calculation.

The Operator Matrix Element
The massive operator matrix element A gg,Q is the expectation value g|O g |g , of the gluonic operator between massless on-shell external gluon states. We will work in R ξ -gauge. Therefore also the corresponding ghost graphs have to be considered. In Eq. (2.1), S and Sp denote the symmetrization of the Lorentz indices and color trace, respectively; F µν is the field strength tensor of QCD and D α denotes the covariant derivative. The OME has been calculated to O(α 2 s ) in [9] and including also terms linear in ε in [10] correcting the previous result.
The renormalized expression of A gg,Q to O(α 3 s ) was derived in [6] and the contributions to O(α 3 s T 2 F N F C F,A ) were calculated in [19]. The OME A gg,Q obeys the expansion with a s (µ 2 ) = α s (µ 2 )/(4π). In the MS scheme with the heavy quark mass m on-shell 3 it is given by gg,Q (2β 0 + 3β 0,Q ) Here δm are expansion coefficients of the renormalization constants for the mass, β i , β i,Q are coefficients of the β-functions (including mass effects), ζ k is the Riemann-ζ function with k ∈ N\{0, 1}, a (2) ij , a (2) ij are two loop contributions to order ε 0 and ε 1 , respectively, and γ ij ,γ ij are the anomalous dimensions. Quantities with a hat in Eq. (2.3) are defined bŷ see Ref. [6]. The unrenormalized OMEÂ (3) gg,Q also receives contributions from the vacuum polarization insertions on the external lineŝ ij,Q are known [5,[8][9][10]13,38]. In particular, all the logarithmic contributions have already been obtained for general values of the Mellin variable N [16,39].
In the following we calculate the O(a 3 s T 2 F C F,A )-contributions to the massive gluonic OME. Before presenting the results, we give a detailed outline of the calculation methods used.

The Methods of Calculation
The gg,Q are given by Feynman graphs with external on-shell gluons (ghosts), a local operator insertion on gluon lines and vertices, and two closed massive quark lines of the same mass m. A calculation along the lines of Refs. [18,20] leads to infinite series which diverge polynomially with degree N. The way to cure this issue will be to separate the variable N from the infinite series by leaving one integral unintegrated. This last integral will then be solved after summation in the space of cyclotomic harmonic polylogarithms [33]. Most of the graphs have been calculated in this way. For a few graphs, we have applied integration by parts and differential equations, see Section 3.5. Throughout the calculation, the results at general values of N are mutually compared to the corresponding moments calculated using MATAD [40].

Feynman Parameterization
The list of graphs was generated with QGRAF [41] and written as momentum integrals using the Feynman rules of [6,42] 4 . The color-algebra was performed using the code Color [44]. The momenta were integrated at the cost of introducing a Feynman parameterization, treating each independent loop separately and introducing for each one of them a family of Feynman parameters. This makes each diagram a linear combination of integrals of the form where for each Feynman parameter family f we used the short-hand notation and ν i are integers denoting the propagator powers. The exponents α i , β i , γ are of the form (a + bε/2) with a, b ∈ Z, and N is the Mellin variable. The operator polynomial is not strictly a polynomial, but in all following cases the δ-distributions and Heaviside functions being present in addition can be removed in such a way that the misnomer is corrected, and the operator polynomial is indeed a polynomial of maximum degree N ∈ N.
The δ-distributions can be integrated using the relations where Y is either a sum of Feynman parameters or a single one. The Heaviside θ-function is defined as These relations are applied in such a way as to keep the operator polynomial as simple as possible. It is indeed possible in all following cases, to map the operator polynomial into one single Feynman parameter, if one uses the following trick: In some cases it is useful to reconstruct a δ-distribution by 4 For the scalar Feynman rules used for the calculation of scalar prototype graphs, see [43].
where X represents a sum of Feynman parameters. Of course the order for the elimination of the Feynman parameters from the θ-functions has to be chosen such that the left hand side of the above equation matches. In this way, an argument (1 − X) consisting of several Feynman parameters is exchanged for only one Feynman parameter. This trick is equivalent to a set of coordinate transformations mentioned in [45] and also used in the calculation of the 2-loop OMEs in [10,13,18,19,46,47]. The above trick has the advantage of giving a clear guideline for how to simplify the polynomial in the N-bracket of the Feynman integrals under consideration. It is worth noting that there are two Feynman parameters, which only occur in the monomial prefactors of the integrand as well as in the operator polynomial. These are due to the fact that the incoming and outgoing momenta are massless. The integral over these Feynman parameters can thus be performed easily, giving simpler N-brackets.
The above methods are applied in order to avoid the proliferation of N. In fact, in all diagrams one can achieve that N only occurs in the exponent of one of the Feynman parameters, allowing to effectively decouple N from the solution of infinite sums. This property of the calculation is of crucial importance, and also carries over to the case of two lines of unequal masses which, however, will be the subject of a future publication.

Mellin-Barnes Representation
The remaining parameters still occur in the denominator polynomial. It has the form (A + B) where A and B are products of elements x i or (1 − x i ), for Feynman parameters x i . Only in the cases of graphs with a massive line that runs through four edges of the graph, e.g. graphs 4 and 5 in Appendix A, a factor (1 − x(1 − y)) in either A or B occurs. A Mellin-Barnes (MB) integral [48,49] is then introduced by the substitution, see e.g. [50,51], This procedure is equivalent to splitting the mass-term off the propagator-like part that occurs in the Feynman parameter representation of a massive vacuum polarization diagram, before proceeding with successive parameterization and momentum integration. In the cases that the products A, B from above factorize completely, all integrals can be performed in terms of Euler's Beta-functions. In the remaining two cases, in which a factor (1 − x(1 − y)) remains, the integrals represent a generalized hypergeometric function 3 F 2 [52,53], which in the scalar diagrams is already given in a form such that it reduces to a ratio of Γfunctions. In the corresponding physical cases, these functions lead to double sums, which can be constructed such that they converge, still keeping N separated from the sums in the way described above.
At this point in all the diagrams only one Beta-function remains that contains both N and ξ. This function is rewritten in terms of a Feynman parameter integral, i.e. for corresponding α and β The reason is that the contour of the Mellin-Barnes integral cannot be closed to a single side. One can see this from two representations. On the one hand, the Beta-function which contains N and ξ has the form , (3.9) so that the denominator drops out of the MB-integral. Hence, if the contour is closed to one side and written as the sum of residues, then, due to its convergence condition, for every set of values of the propagator powers there is an N 0 so that for N > N 0 the sum is divergent.
On the other hand, if the Beta-function is written as a Feynman parameter integral over x, then the factor occurs in the integrand. Here a distinction is necessary between values x < 1 2 for which the contour may be closed towards ξ → ∞, and values x > 1 2 for which ξ → −∞ is the convergent choice. For simplicity, we change the order of the ξ-integration such that the contour can be closed to the right in all cases.
After that the quantity raised to the power ξ is mapped onto a single integration variable T Now it is obvious that all contours have to be closed to the right before applying the residue theorem.
It is worthwhile having a look onto convergence issues of the procedure described so far. First, the Mellin-Barnes integral is introduced in the integrand of the multiple Feynman parameter integral. Employing the nomenclature of [51], the contour follows the usual requirement that left-poles (poles of functions Γ(· · · − z)) are to the left of the contour, and right-poles (poles of Γ(· · · + z)) are to the right of the contour. If left-and right-poles are interleaved on the real axis, the contour winds around them separating the two types of poles.
Of course, contours of the above kind can only be found, if the right-poles are separated from left-poles. In cases where this is not obviously the case, we enforce such a separation by introducing a regularization parameter in a consistent manner throughout the Feynman diagram. So it is most convenient to keep symbolic propagator powers from the beginning, and to use substitutions of these symbolic quantities for the introduction of regulators. We will see later at which point the expansion into a Laurent series in these parameters can be performed most conveniently.
The classical procedure for calculating Mellin-Barnes integrals in particle physics proceeds by deforming the contour and subtracting a finite number of residues, such that the remaining contour integral represents a regular function in ε [51,[54][55][56]. In that case, the expansion can be performed on the integrand level, which simplifies the integrand such that Barnes lemmas are applicable. However, since factors of T ξ occur in the arguments of the contour integrals, cf. Eqs. (3.8, 3.11), no Barnes lemmas [49] can be applied 5 .
In the present calculation, it appears more suitable to write down the sums of residues and generate the necessary simplifications and algebraic relations by symbolic summation methods implemented in the package Sigma [31], equipped with suitable limit procedures for infinite sums.
When residues are calculated and the corresponding sums are written down, one has to perform a Laurent expansion in the regularization parameters. Here it is important to observe the singularity structure.
One therefore brings the Γ-function arguments to a standard form, such that all of them are positive for vanishing regulators Here x and ⌊x⌋ represent the fractional and integer parts of the variable x, respectively. The regulators are assumed to be small enough, such that they only contribute to the fractional part. The Heaviside functions are removed by commuting them with summation operators. This can be done using the following operator relations Once the θ-functions are free of any summation parameters, they can be evaluated. Note that they are also free of the Mellin variable N, since it had been separated from the sums by construction.
Once the Γ-functions have been reflected such that the integer parts of their arguments are positive using the relation [53] their expansion in the artificial regulators is straightforward. Yet one additional preparation is necessary for the expansion in the dimensional regulator ε, since the Feynman parameter integrals may not be well defined in the Lebesgue sense for 0 < ε < 1, but rather as an analytic continuation in ε → 0. The expressions are of the form which only converges if ε > a − 1. Nevertheless, using integration by parts, one can shift this integrand such that it is integrable for 0 < ε < 1. For the form above with a ≥ 1, the relation has to be iterated (a − 1)-times. Here the function g(x) must have sufficiently many regular derivatives on [0, 1], which is indeed the case for the integrals in question. Then the integral represents a regular function in ε, the integrand is measurable for 0 ≤ ε < 1 and thus the Taylor expansion commutes with the integration.
Finally, the expansion of the sums in the dimensional regulator ε can be done using the package EvaluateMultiSums. It also manages the call of Sigma routines and performs limits of many expressions. In particular, the package SumProduction was used to condense the huge expressions into a tractable number of compact but larege sums and automatically apply the summation technologies to obtain the final results.
The result of expansion and summation yields an expression which still depends on one integration variable T , and which contains S-sums [34,58] and cyclotomic S-sums [33] of this variable. They can be converted into (cyclotomic) harmonic polylogarithms (HPL) [33], e.g.
The conversions to iterated integrals are performed using ideas of Ref. [33] and applying automated routines of the package HarmonicSums [32][33][34].
The conversion returns iterated integrals evaluated at 1, but with letters that depend on the remaining integration variable. They will be denoted by where f α is a letter from a cyclotomic alphabet Therefore a procedure is needed that maps the class of iterated integrals appearing here onto (cyclotomic) HPLs with the integration variable in the argument. There is such a procedure which was used for deriving properties of two-dimensional HPLs [59] and in the method of hyperlogarithms [20,21,60]. It makes use of the fact that differentiation of a certain type of iterated integrals with respect to variables appearing in the index leads to a drop in the weight of the function, e.g.
In this way the problem can be traced back to properties of rational functions, solving the problem recursively at a lower weight and integrating again over x, where in each recursive call a constant has to be determined. However, in the case of letters containing polynomials of degree 2 or more, this procedure is not applicable directly, since in general the weight does not drop due to differentiation, e.g. .
Here the following procedure will be useful. Let us distinguish the letters using indices α and denote the corresponding rational functions with f α (x). One can form new letters by scaling the argument of the rational functions with a variable y, cf. Eq. (3.18). If one such letter is built into a cyclotomic HPL with argument x = 1, there is an algorithm for removing the parameter y from the index, such that it occurs in the argument.
At first, by virtue of the shuffle algebra, the weighted letter is brought to the right-most position. Then indexing general rational letters with α i , i = 1, ..., n, we find the algorithm where a partial fraction decomposition is performed in the last step. After that the formula may be recursively applied where the final step is obviously So the result is a multivariate polynomial in iterated integrals of arguments 1, y. This procedure produces also letters 1/(1 − x), which introduce branch points at y = 1. However, considering the integration contour of the iterated integrals infinitesimally away from the real axis does not affect the algorithm introduced above. In this sense, the iterated integrals may be analytically continued, as described in [61] and implemented in HarmonicSums [32][33][34]. Thus they can always be expressed as iterated integrals with arguments in [0, 1].

The Final Integral
Once the sums are performed, i.e. written in terms of iterated integrals, the remaining task is to perform the last integration, which carries the nontrivial dependence on N. However, the integral does not just represent a Mellin-transform, but it contains rational functions R(T ) ∈ {1/(1 + T 2 ), T 2 /(1 + T 2 )}, which are raised to the power N. It therefore seems most natural to consider the generating function of the sequence in N, and to introduce the corresponding tracing parameter κ in the following way .
The integral over T from 0 to 1 is performed in two steps : First a primitive is calculated for the integral in terms of iterated integrals. Then the limits T → 1 and T → 0 are computed. This procedure introduces additional letters into the otherwise cyclotomic alphabet of HPLs, namely Obviously this leads again to re-scaled letters, and one can use the algorithm from above to transform the emerging cyclotomic HPLs at 1 with weighted letters into cyclotomic HPLs with unweighted letters and a function of κ in the argument. It is not hard to see that the functions occurring in the arguments of these HPLs are the functions g(κ) from above.
The limit T → 0 has to be taken carefully to cancel factors of 1/T . Therefore a Taylor expansion is performed. In many cases, relations similar to Eq. (3.17) are used, reading them from right to left in order to obtain the Taylor series. However, such relations are not implemented in HarmonicSums for the additional (weighted) letters of Eq. (3.25). This is due to the requirement of special assumptions on the values of g(κ). We rather use an easy trick to obtain the Taylor series of cyclotomic HPLs extended by the above letter, using the fact that the above letter can be factorized over the complex numbers Then these linear letters are treated like the letter from the alphabet of multiple polylogarithms [34], treating a as real and positive. For the cyclotomic HPLs extended by one such letter, the Taylor series expansions can be derived [33,34]. Finally, the imaginary factors i g(κ) are re-substituted. The results are checked to be regular at T = 0 and thus the limit can be taken. Once the last Feynman parameter integral is performed, we need to find the N th coefficient of the Taylor expansion in κ. For this we would like to make use of methods applicable to HPLs and cyclotomic HPLs which are implemented in the package HarmonicSums [32][33][34]. It is therefore necessary to make sure that the dependence on ln(κ) cancels. These terms can be eliminated using argument transformations and algebraic relations [62] of the (cyclotomic) HPLs. 6 At first, the arguments are mapped back into the interval [0, 1] where the length of the list β is bounded by the length of α, and the a β are integer coefficients.
Relations of this kind can be obtained algorithmically and are implemented for all cyclotomic HPLs in the package HarmonicSums.
Then the square roots are removed from the arguments, as far as possible. For this step one makes use of the fact that all cyclotomic HPLs with arguments x 2 can be rewritten in terms of cyclotomic HPLs with arguments x. These transformations can be inverted, so that (cyclotomic) HPLs which contain the letter f (1,0) (x) = 1 x−1 and the argument √ 1 − κ are mapped onto (cyclotomic) HPLs with argument 1 − κ and (cyclotomic) HPLs without the letter f (1,0) , i.e.
where in the vector α there is an index (1, 0). The length of β is again bounded by the length of α, and γ is free of the index (1, 0). This reduction is, however, not complete so it is introduced by constructing a basis of HPLs w.r.t. the shuffle relations as well as the relations of squared arguments. It is a sign of a proper Laurent-series that after the reduction to such a basis the remaining (cyclotomic) HPLs involving the letter f (1,0) and with argument √ 1 − κ will cancel. The last step to properly cancel logarithmic parts is to write all ln(κ) parts explicitly, using the flip relation In the present case, this relation has to be applied only to HPLs with letters from the alphabet This subset is not closed under the flip x → (1 − x), so the property will lead to multiple polylogarithms [34] in the result. Nevertheless, the representation is standardized so that indeed all dependencies on ln(κ) cancel. The remaining HPLs fit into the alphabet where the letters f 0 , f −1 , f (4,0) , f (4,1) occur in HPLs with arguments √ 1 − κ, and letters f 1 , f −1 , f (4,0) , f (4,1) lead to HPLs with argument κ. 7 The result thus obtained has a Taylor expansion in κ around 0. The remaining step to obtain the all-N result is to extract the N th coefficient of the corresponding Taylor series. This can be done analytically term by term, using expansions of individual factors and calculating their Cauchy products, as well as by deriving difference equations which are solved in terms of indefinite nested sums. Also these methods are available through the packages HarmonicSums and Sigma.
As a result of this procedure one obtains a large expression in terms of sums of higher depth, involving definite and indefinite sums and products. To obtain a minimal representation, the package EvaluateMultiSums and Sigma can be applied, in order to represent these objects in terms of indefinite nested sums, and in order to eliminate all relations among these indefinite nested sums and products to obtain a basis-representation.

Operator Insertions on External Vertices
The class of graphs with two massive fermion lines of the same mass also includes graphs with operator insertions on external gluon vertices. In the scalar case these graphs are directly related to graphs with operator insertions on a line, see [20] for similar properties used in the calculation of ladder graphs.
The idea carries over to the physical case, but there are no simple relations among graphs. Instead if a method is known for the calculation of certain graphs with operator insertions on lines, then the same methods apply for the graphs with operator insertions on external gluon vertices.
The reason lies in the structure of the Feynman rule for the operator insertion of a gluon vertex, which can be taken from [6,42] V abc µνλ (q 1 , q 2 , q 3 ) = − ig In this notation, the summands in the left column of Eq. (3.34) all behave like operator insertions on lines. Furthermore, if q 1 = p is the external momentum then the first and last summand in the second column behave like insertions on lines too, but here in addition the result is subject to a finite sum of the form The remaining summand (second term, right column) can be summed on the level of Feynman rules, and using q 2 + q 3 = −q 1 = −p one finds In this way, the operator insertion on an external vertex is related to operator insertions on internal lines. However, a direct relation between a graph with a vertex insertion and the corresponding graphs with line insertions does not follow from this consideration, due to the presence of the tensors t 3g µνλ and τ 3g µνλ .

Integration by parts and differential equations
The diagrams shown in Figure 1 turned out to be too cumbersome to be calculated with the methods described before. For this reason, these diagrams were computed using a different approach. For each diagram, a Form program [64] was written in order to replace the propagators and vertices from the output of QGRAF [41] by the corresponding Feynman rules. Further it introduces the corresponding projector for the Green's function under consideration and performs the Dirac-algebra in the numerator. After this, each diagram ends up being expressed as a linear combination of scalar integrals, which were then reduced using integration by parts to master integrals using the program Reduze2 [65] 8 . This is a C++ program based on Laporta's algorithm [68], and has been adapted to the case were we have operator insertions in the integrals. Figure 1: Diagrams calculated using differential (difference) equations. Here the masses for both fermion loops are equal.
In total, sixteen master integrals were needed in order to calculate these diagrams. Eleven of them have the general form where we use the shorthand notation and The superscript D has been included in J D ν 1 ,...,ν 9 (N) in order to make explicit the dependence on the dimension D. The eleven integrals of this type are then The other five master integrals are , (3.54) , (3.55) . (3.56) Notice that integrals J D 15 and J D 16 are just constants w.r.t. N. The integrals J D 12 , J D 13 and J D 14 yield Feynman parameter integrals that can be performed in terms of Beta-functions. We obtain where we have set the mass m, ∆.p and spherical factors to 1 for simplicity. Any given scalar integral will be written as a linear combination of these master integrals. Since the coefficients of these linear combinations may contain poles in ε = D − 4, the master integrals may need to be expanded to higher orders in ε accordingly, in order to get the corresponding scalar integrals up to order ε 0 .
The integrals J D 1 (N), . . . , J D 11 (N) were calculated using the differential equations method [69]. This method has been applied successfully to many problems where Feynman integrals depending on one or more invariants appear. The idea is to take derivatives of the master integrals with respect to these invariants and re-express the result in terms of the master integrals themselves. This leads to a system of differential equations that can then be solved. In the present case, the integrals depend on the invariants m 2 and ∆.p. However, the dependence of the integrals on these invariants is trivial. They are just proportional to (∆.p) N and (m 2 ) −ν+ 3 2 D . Here ν is the sum of powers of propagators. Therefore, taking derivatives with respect to these invariants does not lead to any new information. The integrals have the form and it is actually the calculation of the function F (N) expanded in ε = D − 4 that is non-trivial.
One might think about taking derivatives with respect to N, but this changes the structure of the integrals in a way that does not allow the application of the differential equations method. In view of this, we introduce a new parameter x and rewrite the operator insertion in the following way (3.60) By doing this, we trade the dependence of the integrals on N by a dependence on x, and the operator insertion becomes a denominator that can be treated as an additional artificial propagator. In fact, it is in this x-representation of the integrals in which all the reductions to master integrals are performed using Reduze2. Laporta's algorithm requires integrals to have definite powers of propagators, and although it may be possible to express ∆.k 3 in terms of inverse powers of propagators (by taking ∆ as an external momentum), one faces the problem that the power N in the operator insertion (∆.k 3 ) N is arbitrary. By turning the operator insertion into an artificial propagator as in Eq. (3.60), we circumvent this difficulty. Let us defineĴ .

(3.61)
We can now take derivatives with respect to x. 9 This will raise the power of the artificial propagator, leading to integrals that can then be reduced, and as usual, a system of differential equations is generated. For example, the first three integrals in Eqs.
where we have set m 2 and ∆.p to 1 for simplicity. This system can now be solved, provided the constant integrals J D 15 and J D 16 are previously computed, and a few initial conditions are provided. These initial conditions will be the values of the integrals and some of their derivatives at x = 0. Since the N th derivative of a given integral J D i (x) at x = 0 is equal to N!J D i (N), we see that giving these initial conditions is equivalent to giving a few initial values for J D i (N). When we take the derivatives of the remaining integrals in Eqs. (3.44)-(3.51), the integrals J D 1 , J D 2 and J D 3 will also appear on the right hand side of the equations. For example . Likewise, J D 4 (x) will appear on the right hand side of the differential equations of the next integrals, etc. We can see that we must solve the system of differential equations starting with the simplest integrals, and gradually incorporate the results to solve the more complicated ones. This is all done with the help of the Mathematica packages Sigma [31], HarmonicSums [32][33][34], EvaluateMultiSums, SumProduction [35], and OreSysG [70]. These packages construct a system of difference equations from the differential equations, and then solve for J D i (N) directly, instead of J D i (x). For example, one may transform the system (3.62)-(3.64) using Eq.
Here we left out the explicit dependence on the dimension D = 4 + ε in the functions J i . Let us discuss now the calculation of the initial values required in order to solve the differential (difference) equations discussed above. These are basically the values of the integrals for a few fixed values of the Mellin variable N. In some cases, these values are needed only up to order ε 0 , which can therefore be obtained using MATAD [40]. More often, the initial values are needed up to higher orders in ε, and a different method to obtain them must be used. In the following, we describe the method we used in such cases based on the α-parameterization of the integrals. In the present calculation, five initial values starting from N = 1 were needed up to order ε 2 for the master integral and two initial values up to order ε starting from N = 1 were needed for J D 7 = J D 2,1,0,1,0,1,1,0,0 (N) = dk (3. 70) In what follows, the masses appearing in some of the propagators do not play any role, so we will omit them for the time being. Let us consider the general integral in Eq. (3.38). Removing the operator insertion (i.e., taking N = 0) the α representation of this integral is given by and q 1 = −β 2 p, q 2 = −β 4 p and q 3 = −β 9 p , (3.73) where the β i 's are defined using the θ-function as The product of integrals in the α parameters, and the sum in the exponential in the first and second lines of Eq. (3.71), run over the values of l, i and j corresponding to the propagators that are actually present in the integral under consideration. We can now introduce the operator insertion in our integrals in the following way, cf. also [51], , (3.75) and now q ′ 1 = −β 2 p, q ′ 2 = −β 4 p and q ′ 3 = −β 9 p + r∆ . (3.76) In the case of integral J D 1 (N), we get and q ′ 1 = −α 2 p, q ′ 2 = 0 and q ′ 3 = r∆ , (3.78) and one has α 6 α 7 α 3 α 6 + α 7 α 6 α 6 α 7 α 2 α 6 + α 7 α 6 + α 2 α 7 α 2 α 7 + α 6 α 7 α 3 α 6 + α 7 α 6 α 2 α 7 + α 6 α 7 If we apply Eq. (3.75) in this case, we obtain Here the operator i + shifts the power of the i th propagator by one, and also multiplies the integral by −ν i , i.e. i + J D ν 1 ,...,ν i ,...,ν 9 = −ν i J D ν 1 ,...,ν i +1,...,ν 9 . The integral J D 16 is more complicated. After Feynman parameterization we obtain (3.92) We can now obtain a Mellin-Barnes representation for this integral by splitting the denominator in the equation above using which leads to In this representation, the integral can be calculated with the help of the Mathematica package MB [56]. This package finds a value for γ and ε = D − 4 such that the integral in Eq. (3.94) is well defined. Then it performs an analytic continuation to ε → 0 and expands in ε. After this, we can close the contour to the right or to the left and take residues. This leads to sums that can be performed with the package Sigma. For the different shifts in D, we obtain  gg,Q of the unrenormalized OME (2.8). Defining it is given by with the polynomials P i (4.6) (4.7) P 6 = 99N 6 + 297N 5 + 631N 4 + 767N 3 + 1118N 2 + 784N + 168, (4.8) 14) Here we use the short-hand notation S a (N) ≡ S a for the harmonic sums [27]. The polynomial denominators in Eq. (4.2) show evanescent poles at N = 1/2, 3/2. However, the function is continuous at these points, as the expansion around these values shows. This also confirms that the rightmost pole is located at N = 1 as expected for a gluonic quantity in QCD. This also applies to the OME, Eq. (4.26). In Eq. (4.2) the new sum  [16] the OME behaves the same way, cf. [73].
It is an interesting question, as to whether new structures, like those in Eq. (4.16) compared to the usual harmonic sums, can be recognized in studying the minimal difference equations 10 they obey. For this purpose we consider the equation for the harmonic sum S 2,1 (N) at one side and Eq. (4.16) on the other side. The former obeys the difference equation The term T without the ζ 3 -contribution obeys with the initial values . (4.21) Both difference equations are of degree and order three and are of quite similar structure. The different type of the solutions are therefore hardly recognized ab initio. The Mellin inversion of the binomial terms yield [28] N j=1  Therefore the two new letters [28] f w 1 (x) = 1 appear in the x-space representation beyond those forming the usual harmonic polylogarithms [61].

The Operator Matrix Element
The O(T 2 F C F,A ) contribution to the operator matrix A gg,Q is given by using Eqs. (2.3, 4.2) and the corresponding expressions implied by renormalization from Ref. [16]. The polynomials P i read (4.33)   The analytic continuation of the OME Eq. (4.26) from the even moments N = 2n, n ∈ N to the complex plane is obtained using the asymptotic representation for the harmonic sums [73,74] and Eq. (4.16) supplemented by the recursion relations for N → (N − 1) of Eq. (4.26). The OME in the MS scheme is obtained by the following transformation with the polynomials (4.47) Here we have set the masses in both schemes equal symbolically, to obtain a more compact expression.

Anomalous Dimension
As a by-product of the calculation we obtain the corresponding contributions to the anomalous dimensions from the single pole term 1/ε or the corresponding linear logarithmic term, cf. Eq. (2.3), where Eq. (4.48) confirms previous results in [38] by a first direct diagrammatic calculation, here for massive graphs containing two fermion lines of equal mass. In Ref. [19] the anomalous dimension has been confirmed for 3-loop graphs containing one massless and a massive fermion line.

Conclusions
The contribution of O(T 2 F C F,A ) to the massive operator matrix element A gg,Q (N) at 3-loop order has been calculated. It receives contributions from diagrams with two internal massive quark lines of equal mass. The OME can be expressed in terms of harmonic sums, supplemented by a single new binomially weighted harmonic sum. The analytic continuation to N ∈ C is given by the recurrence relation of the expressions and the asymptotic representation. The OME has poles for N ∈ Z, N ≤ 1. The results have been given for both the on-shell and MS-scheme for the heavy quark mass. In the latter scheme, terms ∝ ζ 2 are not present, cf. also Ref. [6]. As a by-product the corresponding contribution to the 3-loop anomalous dimension γ gg has been obtained in an independent calculation ab initio. The calculation of the diagrams with two massive fermion lines need more special techniques than in the case of a single fermion line. Here the use of Mellin-Barnes representations and generating functions based on cyclotomic harmonic polylogarithms and S-sums is essential. In some of the diagrams we applied the method the integration by parts method and applied differential equations to calculate the associated master integrals. The technologies described can be generalized to the case of two different masses.

A Results for the scalar graphs
In the following, the results for the scalar prototypes of the graphs contributing to the O(T 2 F C F,A ) part of the operator matrix element A (3) gg,Q are summarized. These diagrams are much simpler to calculate than the corresponding complete diagrams. However, they show the principal structures of the full diagrams and share a common calculational scheme. The large amount of numerator terms and their variation, however, increases the complexity of the QCD diagrams significantly.
All diagrams are normalized such that the factor In some of the graphs, the denominator structure show evanescent poles at N = 1/2, 3/2, 5/2. The expansion of the whole function around these values shows continuity.