The 3-Loop Pure Singlet Heavy Flavor Contributions to the Structure Function $F_2(x,Q^2)$ and the Anomalous Dimension

The pure singlet asymptotic heavy flavor corrections to 3-loop order for the deep-inelastic scattering structure function $F_2(x,Q^2)$ and the corresponding transition matrix element $A_{Qq}^{(3), \sf PS}$ in the variable flavor number scheme are computed. In Mellin-$N$ space these inclusive quantities depend on generalized harmonic sums. We also recalculate the complete 3-loop pure singlet anomalous dimension for the first time. Numerical results for the Wilson coefficients, the operator matrix element and the contribution to the structure function $F_2(x,Q^2)$ are presented.


Introduction
The present precision of deep-inelastic scattering data [1] allows for the measurement of the strong coupling constant α s (M 2 Z ) at an accuracy of 1% [2] and for a precision determination of the parton distribution functions [3,4] and the mass of the charm quark m c [5]. In the future, new dedicated deep-inelastic experiments may be carried out at high luminosities at the EIC [6] and at even higher energies than those at HERA [7] as planned for the LHeC [8]. At these facilities the experimental resolution will be even higher. The corresponding analyses require 3-loop accuracy, including the heavy flavor Wilson coefficients. At present the heavy flavor corrections to deep-inelastic scattering are known to 2-loop order in semi-analytic form [9] 1 . As has been shown in Ref. [11] for scales Q 2 > ∼ 10 m 2 , with m the heavy quark mass, the heavy flavor contributions to the structure function F 2 (x, Q 2 ) can be calculated to the 1% level employing a factorization of the scattering cross section into massive operator matrix elements (OMEs) and massless Wilson coefficients [11]. This enables us to calculate the higher order corrections in Quantum Chromodynamics (QCD) in analytic form.
At 3-loop order this calculation has been performed for a series of Mellin moments in Ref. [12] in 2009. The calculation of the corresponding results for general values of the Mellin variable N requires by far different techniques than those having been used in [12]. To a wide extent, they were not previously available and had to be newly developed in course of the present calculation. In the past we have recalculated and corrected the 2-loop results [11,[13][14][15] using more systematic summation and integration methods in Refs. [16][17][18][19][20]. Furthermore, we calculated the asymptotic heavy flavor corrections to the structure function F L (x, Q 2 ) [21,22]. Very recently, we presented results on the operator matrix element A (3) gq [23] and the flavor non-singlet OMEs and Wilson coefficients [24]. Furthermore, the 3-loop contributions of O(N F T 2 F ) have been computed completely [25,26], as well as the contributions O(T 2 F ) to the OMEs A gq and A gg [27] stemming from graphs with two internal fermion lines carrying the same mass 2 . Technical aspects of these calculations have been presented in Refs. [29,30]. In all these calculations the respective contributions to the 3-loop anomalous dimension are obtained as a by-product.
In the present paper we calculate the pure singlet contributions to the heavy flavor Wilson coefficient H PS 2,q at 3-loop order in the asymptotic region and present the operator matrix element A , which also appears as one of the matching coefficients in the variable flavor number scheme (VFNS). As in previous calculations [27,31], new mathematical structures emerge in intermediary steps. In x-space they appear as generalized harmonic polylogarithms [32]. In the physical result they can be mapped back to the usual harmonic polylogarithms [33] at the arguments x and a new one at y = 1−2x. In Mellin-N space, generalized harmonic sums contribute to the result [34].
The paper is organized as follows. In Section 2 we briefly describe the basic formalism. Some technical aspects of the calculation of the massive OME are outlined in Section 3. This concerns the reduction of the Feynman diagrams to master integrals and the different methods we have applied for their calculation. The unpolarized pure singlet anomalous dimensions to 3-loop order is presented in Section 4 and compared with results in the literature. The massive OME A is presented in Section 6 and numerical illustrations are given for the pure singlet contribution to the structure function F 2 (x, Q 2 ) due to charm and bottom quarks. Section 7 contains the conclusions. In Appendix A we discuss aspects of the contributing integral families. Mellin representations of the newly contributing generalized harmonic sums are given in Appendix B. The expression for the operator matrix element A PS Qq and the asymptotic massive Wilson coefficient H PS 2,q in x-space are given in Appendix C.

Basic Formalism
The renormalized pure singlet OME in the MS-scheme for the coupling constant to 3-loop order [12] is given by It describes the transition between massless on-shell quark states q|, characterized by a local quark operator in the light-cone expansion [35], which is located on the heavy quark line. The corresponding pure singlet contributions in case the operator is located on an internal light quark line has been dealt with in Refs. [22,25]. The OMEs at 2-and 3-loop order are given by (2),PS Qq (β 0 + β 0,Q ) + 8γ (0) qg a (2) gq,Q − 8γ (0) gq a (2) Qg −γ (0) qg γ (0) gq ζ 2 γ (0) gg − γ (0) qq +6β 0 + 8β 0,Q ln m 2 µ 2 + 4(β 0 + β 0,Q )a (2),PS Qq ij , k = 1, 2, 3 is the constant part of the unrenormalized OME at O(a k s ), with the strong coupling constant g s expressed as a s = g 2 s /(4π) 2 ≡ α s /(4π),ā (k) ij , k = 1, 2 denotes the part ∝ ε of the unrenormalized OME at O(a k s ), with ε = D − 4 the dimensional parameter, β k and β Q,k are the expansion coefficients of the QCD β-function in the MS-scheme and for massive contributions, δm (l) k are the expansion coefficients of the renormalized quark mass m, µ is the renormalization scale, N F denotes the number of light quark flavors, and ζ k = ∞ l=1 (1/l k ), k ∈ N, k ≥ 2 denotes the Riemann ζfunction at integer argument. For details of the notation see Ref. [12]. Here and in the following we also use the shorthand notationŝ In the asymptotic region Q 2 ≫ m 2 the pure singlet heavy flavor Wilson coefficient is given by [12] H Here C (l) 2,j , with j = q, g, l = 1, 2, 3 denote the corresponding light flavor Wilson coefficients [36][37][38][39]. The OME A (2),PS Qq has been calculated in [11,16] and A (2) gq,Q in [14,18].
The heavy flavor pure singlet contribution to the structure function F 2 (x, Q 2 ) in case of the coupling of the exchanged virtual photon of virtuality Q 2 to the heavy quark line of charge e Q is obtained by [22] F QQ,PS 2 (x, Q 2 ) = e 2 Q x 1 x dy y H PS 2,q y, The quark-singlet distribution is given by and q(x)(q(x)) denote the light flavor quark and anti-quark number densities, respectively. Before we present the physical results on the pure singlet 3-loop anomalous dimension, the massive OME and Wilson coefficient, and numerical results on the pure singlet contribution to the structure function F 2 (x, Q 2 ), we discuss a series of technical details of the present calculation.

Details of the Calculation
The massive OME A (3),PS Qq is represented by 125 Feynman diagrams, a sample of which is shown in Figure 1. The diagrams are generated using QGRAF [40]. Here the operator insertions are realized in terms of vertices with non-propagating scalar particles attached to them, cf. [12]. The propagators, vertices and operator insertions from the output of QGRAF are then replaced by the corresponding Feynman rules using a FORM [41] program [12], which also allows us to introduce the corresponding projector for the Green function under consideration and perform the Diracmatrix algebra in the numerator of the Feynman integrals. After this, the diagrams end up being expressed as linear combinations of scalar integrals. In the case of the contributing bubble topologies we used hypergeometric techniques [16][17][18][19][20][42][43][44] and calculated the corresponding graphs directly, cf. [25,45]. The packages Sigma [46,47], EvaluateMultiSums, SumProduction [48], ρsum [49], HarmonicSums [32,50] and OreSys [51] have been used extensively. There are nine types of 3-loop integrals involved in the calculation of A  Here ∆ denotes a general light-like vector. The Feynman rules, including those for the local operator insertions, are given in Ref. [12]. . The dashed arrow lines represent massless quarks, while the solid arrow lines represent massive quarks, and curly lines are gluons. In terms of Feynman integrals, diagrams (a-j) represent the main topologies, in other words, other diagrams, such as diagrams (k)-(t) can be seen as sub-topologies and/or are related to diagrams (a-j) by symmetry. The symbol ⊗ denotes the local operator insertion, see Ref. [12].
In Eqs. (3.1-3.5), ν 1 , . . . , ν 9 , a, b and c are integers, and we use the shorthand notation The (inverse) propagators D 1 , . . . , D 9 are given by where m is the mass of the heavy quark and p is the momentum of the external massless quark, which is taken on-shell (p 2 = 0). For example, the diagram in Figure 1a can be written as a linear combination of the K 1 -type integrals defined in Eq. (3.1), and the diagram in Figure 1b can be written in terms of the K 3 -type integrals defined in Eq. (3.3).
The last four types of integrals have a different propagator structure and are given by (3.10) For example, the diagrams in Figures 1e, 1f and 1g can be written as linear combinations of the integrals defined in Eqs. (3.8), (3.10) and (3.9), respectively. The Feynman rule for 4-point operator insertions contains two terms. For some diagrams, such as those in Figures 1h and 1i, both parts of the diagram associated with each term can be written as a linear combination of the K 5 -type integrals defined in Eq. (3.5). In the case of the diagram in Figure 1j, one piece of the diagram can be written in terms of the K 5 -type integrals, and the other piece in terms of the K 9 -type integrals defined in Eq. (3.11). Any given diagram has at most eight propagators, so at least one of the propagators in the lists D 1 , . . . , D 9 or D ′ 1 , . . . , D ′ 9 plays the role of an auxiliary propagator, whose presence allows us to uniquely express all possible scalar products of momenta k i · k j and k i · p (i, j = 1, 2, 3) as linear combinations of all inverse propagators D 1 to D 9 (or D ′ 1 to D ′ 9 ). Which one(s) of the nine propagators turn out to be auxiliary depends on the specific diagram under consideration. Scalar integrals will be identified by the indices ν 1 , . . . , ν 9 , a, b and c, where some of the indices ν 1 to ν 9 can be negative, which will represent a scalar integral with irreducible numerators. The factors (∆.k 1 ) a , (∆.k 2 ) b and (∆.k 3 ) c arise from contractions of an internal momentum with a ∆ appearing in the operator insertion Feynman rule. In the case of integrals (3.1), (3.2) and (3.9), the indices a, b and c are bounded by 0 ≤ a + b + c ≤ 1, while in the case of integrals (3.3), (3.4), (3.8) and (3.10), we have that 0 ≤ a + b + c ≤ 2, and in the case of Eqs. (3.5) and (3.11) we get 0 ≤ a + b + c ≤ 3.

Integration by parts identities
The number of scalar integrals required in order to calculate the diagrams is quite large. We use integration by parts identities [52] in order to express all scalar integrals in terms of a much smaller set of master integrals. For this purpose we use Reduze2 [53] 3 , which is a C++ program based on Laporta's algorithm [56][57][58][59]. It is somewhat difficult to adapt this algorithm to the case where we have operator insertions since it requires the integrals to be identified by definite indices, and in the numerator of the integrals we have dot products of internal momenta with ∆ raised to arbitrary parameters such as N , j or N − j. For this reason, we introduce a generating function in a new variable x, rewriting all operator insertions in terms of a sum in N , cf. [29]. For example, a fermion line insertion with momentum k going through the line will be re-expressed as 4 This then can be treated as an additional propagator 5 , and Laporta's algorithm can be applied without further modifications. Similarly, the 3-point and 4-point vertex operator insertions can be replaced by products of two or three such artificial propagators, respectively. In the 3-point case, we get and in the 4-point case, the replacement holds. In this way, the five integrals given in Eqs. (3.1-3.5) can all be represented in terms of the general integral J B1a ν 1 ,...,ν 12 (x) = dk where (3.17) and the propagators D 1 to D 9 are the ones defined in Eq. (3.7). Notice that the set of (inverse) propagators D 1 to D 12 is complete and minimal, which means that any scalar product of a loop momentum with ∆, p or loop momenta can be uniquely expressed as a linear combination of these propagators. A set of propagators satisfying this condition is called an integral family. The superscript B1a in Eq. (3.16) labels the particular integral family defined by the propagators D 1 to D 12 . A given scalar integral will be completely identified by specifying the integral family and the set of indices ν 1 to ν 12 . There is a total of 24 integral families needed for the calculation of all operator matrix elements, although, as we will see, only three of them are needed for the calculation of A K 2 ({ν i }; 0, 0, 0; N ) → J B1a ν 1 ,...,ν 9 ,1,0,0 (x), (3.19) ..,ν 9 ,1,1,0 (x), (3.20) ..,ν 9 ,1,0,1 (x), (3.21)  (3.20), respectively. The diagrams on the lefthand side of Figure 2 must be interpreted as the corresponding scalar integral with no numerator other than the term coming from the operator insertion shown below the diagram. The diagrams on the right-hand side of Figure 2 represent the scalar integrals after the transformations in Eqs. (3.13-3.15) are done. Solid and dashed lines represent massive and massless propagators, respectively. A large dot on a line in these diagrams represents an artificial propagator of the form (1 − x∆.q) −1 , where q is the momentum going through the line in the depicted direction. If any of the indices a, b, c in Eqs. (3.1-3.5) is different from zero, one can always represent the corresponding numerators via the propagators D 10 , D 11 and D 12 , leading to linear combinations of the integrals defined in Eq. (3.16). For example, ..,ν 9 ,0,0,0 (x) − J B1a ν 1 ,...,ν 9 ,1,−1,0 (x) (3.23) or Similarly, the integrals defined in Eqs. (3.8-3.9) can be re-expressed as where 26) and the integrals defined in Eq. (3.10) and (3.11) are associated with We can see that ..,ν 9 ,1,0,1 (x), (3.29) ..,ν 9 ,1,0,0 (x), (3.30) ..,ν 9 ,1,1,0 (x), (3.31) ..,ν 9 ,1,1,1 (x), (3.32) and in cases where any of the indices a, b, c is different from zero, we again get linear combinations similar to those of Eqs. (3.23-3.24). The integral families B1a, B5a and B5c are shown again in Appendix A, where we depict the different topologies that they cover. In total, 66 master integrals were required for the reduction of all integrals appearing in the calculation of A . Of those, 55 belong to family B1a and 11 to family B5a. In Table 1, we list the master integrals in family B5a. The list of integrals in family B1a is a bit long and will be omitted here. No master integrals in family B5c were required, since all integrals in this family were reduced to master integrals in family B5a. This is a peculiarity of A . For other operator matrix elements where family B5c appears, a few master integrals belonging to this family will be required.
It must be pointed out that the basis of master integrals we have chosen is arbitrary, and we can in principle choose any other basis. For convenience, we have chosen a basis where no master integral with negative indices appear. This choice was motivated by the fact that this type of integrals are easier to handle for many of the methods we have used to solve them. These methods are discussed in the next section.

Calculation of the master integrals
For the calculation of the master integrals we used a variety of methods. For the simplest cases, we combined the propagators using Feynman parameters, leading to expressions that can be solved in terms of generalized hypergeometric functions [42][43][44] or by introducing a Mellin-Barnes [60] representation. For more complicated integrals, we used the differential equations method. Below we will describe these methods using a few illustrative examples.
As we have seen, the introduction of the variable x has allowed us to turn all operator insertions into artificial propagators, making the application of Laporta's algorithm straightforward. The calculation of the integrals in this representation using Feynman parameters, although possible in principle, can be somewhat difficult, since integrals become more complicated as more propagators are present (more Feynman parameters need to be introduced). For this reason, in these cases we calculate the master integrals in the original N -dependent representation. Once the integrals are calculated as functions of N , one can always go to the x-representation when needed by performing the transformations given in Eqs. (3.13-3.15). On the other hand, in the case of the differential equations method, we will see that the introduction of the variable x turns out to be actually quite advantageous, although this method ultimately leads to difference equations in the variable N , and we end up obtaining the integrals as functions of N , just like in the other methods. In calculating the master integrals and for their assembly to the individual Feynman diagrams we made also use of the package Matad [61] and have performed checks for fixed moments.

Hypergeometric functions and summation methods
The majority of the master integrals were calculated in terms of hypergeometric functions evaluated at 1, or multiple sums of such functions where the summation indices, the dimensional parameter ε = D − 4 and N may appear in the parameters of the function. If the corresponding series representation is convergent, the resulting multiple sums can then be evaluated with the Mathematica packages Sigma, HarmonicSums, EvaluateMultiSums and SumProduction. These packages implement summation algorithms based on difference fields [62][63][64][65][66][67][68][69][70] and can deal with finite and infinite sums, simplifying the expressions in terms of definite nested sums and products.
Let us consider, for example, the following master integral After we introduce Feynman parameters and perform the loop momentum integrals, we obtain the following expression The integral in w gives just a Beta-function, while the integral in x can be done in terms of a hypergeometric 2 F 1 function. We get, We can now use the following analytic continuation [71], which leads in our case to Now we can do the shift y → 1 − y, and the remaining integrals in y and z can be evaluated using [43] 1 We obtain an expression in terms of 4 F 3 hypergeomteric functions evaluated at 1, The parameters of the hypergeometric functions above satisfy the criteria for convergence, so we can use the corresponding series representations. We find .

(3.41)
So we get an expression in terms of a double sum (one of them finite and the other infinite). This double sum can be done using the packages we mentioned at the beginning of this section. These packages can perform the ε expansion of the expression given above to the required order, and then calculate the sums. The final result is Here S a (N ) denote the nested harmonic sums [72,73]. They are defined by We use the shorthand notation S a (N ) ≡ S a . Note that we have omitted an overall factor of iS 3 ε , where S ε is the spherical factor given by Here γ E denotes the Euler-Mascheroni constant. The expression given in Eq. (3.42) is divergent for N = 1, so we have to calculate this value separately. We get,

Mellin-Barnes integral representations
A few master integrals were calculated using a Mellin-Barnes integral representation. In particular, we used this method for seven K 7 -type master integrals, corresponding in the x-representation to the integrals in Table 1 starting from the third row until the ninth row, together with the first integral appearing in this Table (which is independent of x). Let us consider the case After Feynman parameterization we obtain Now we split the last term using, cf. [60]. This makes it possible to integrate the Feynman parameters at the expense of introducing the contour integral in σ. We obtain, . (3.49) At this point we use the Mathematica package MB [74] to find a value for γ and ε such that the integral in Eq. (3.49) is well defined, and then analytically continue to ε → 0 and later expand this expression in ε. We get, where a 0 (N ) is a term produced by MB after taking a residue at σ = 3 2 ε in order to perform the analytic continuation in ε. It is given by The functions b 0 (N ) and b 1 (N ) are the following contour integrals, Here ψ(z) denotes the Digamma function. We can now close the contours to the left (or to the right), express the integrals in terms of a sum of residues and obtain The above sums can now be performed using the Mathematica package Sigma. The final result is The generalized harmonic sums [32,34] are defined by and we use the shorthand notation S a ( b, N ) ≡ S a ( b). Here we have again omitted an overall factor of iS 3 ε defined in Eq. (3.44). The expression in Eq. (3.56) is divergent for N = 1, so this value has to be computed separately. We obtain The constant B 4 is given by with σ a = lim N →∞ S a (N ), and belongs to the multiple zeta values [75]. Here Li n (x) denotes the polylogarithm [76].

Differential equations
In the x-representation of the integrals, we have the possibility to take derivatives of the integrals with respect to x. If we do this to a master integral J(x, ε), the result can then be rewritten using integration by parts (IBP) reductions in terms of the master integrals J i (x, ε) themselves.
(i.e., integrals for which the set of indices ν i that are positive is the same) will produce a system of coupled differential equations, which we can solve after an expansion in ε. In Table 2, we show the list of integrals solved using this method. They all belong to the integral family B1a. We have included a few horizontal lines separating the different sectors. Let us consider the first two integrals in Table 2, and use the following shorthand notation, Taking derivatives with respect to x we obtain where T 1 (x) and T 2 (x) are linear combinations of sub-sector master integrals that have been solved previously. T 1 (x) and T 2 (x) can be turned into the N -representation using Sigma. Then using the fact that we get the following system of coupled difference equations: whereT 1 (N ) andT 2 (N ) are the N th terms of the Taylor expansions of T 1 (x) and T 2 (x), respectively. The system can be solved using the Mathematica packages Sigma [46,47] and OreSys [51]. In order to do so, we need to provide a few initial values. This can be done along the lines of Ref. [27]. We obtain The results are given in terms of standard harmonic sums. We give the results up to ε 0 , although F 1 (N ) and F 2 (N ) are needed and were calculated up to orders ε and ε 2 , respectively. We obtain Here again we have omitted an overall factor of iS 3 ε in both expressions. Further details on our differential equation method are outlined in [77].

The Pure Singlet Anomalous Dimension
The pure singlet anomalous dimension at 3-loop order can be calculated in complete form from the term ∝ ln(m 2 /µ 2 ) of the renormalized OME Eq.
as shorthand notation. In the present calculation we obtained for the 2-and 3-loop anomalous dimensions with the polynomials Here the color factors are given by The three-loop pure singlet anomalous dimension depends on the following harmonic sums only if one reduces the final result algebraically [78], cf. also [79]. These sums, furthermore, obey structural relations [80,81], leading to a further reduction to In x-space the pure singlet anomalous dimensions read : Here we used the shorthand notation H a (x) ≡ H a for the harmonic polylogarithms [33]. They are defined by The letters f b (x) are given by Furthermore, the kth iteration of the letter 0 leads to ln k (x)/k!. Again, we have used the algebraic relations [78]. The pure singlet anomalous dimension up to 3-loop order depends on the following harmonic polylogarithms only. Since the functions H 0,−1,1 and H 0,1,−1 emerge as a sum, all these harmonic polylogarithms can be represented as Nielsen integrals [82][83][84], with argument x, −x and x 2 , respectively, cf. [23,80]. These are the functions This behaviour is observed for all 3-loop splitting functions, see Refs. [80,85]. Despite of the algebraic reduction, the representations in x-space request more basic special functions than the case in Mellin-N space. Eqs. (4.2, 4.19) agree with the results given in Refs. [38,86] and Eqs. (4.3, 4.20) with the corresponding results given in [87]. For the latter case we present the first independent recalculation here.
It is related to the Euler-Mac Laurin [91] representation of linear combinations of harmonic polylogarithms and can be extended to arbitrary precision. It improves earlier representations based on Chebyshev polynomials [92] as applied to Nielsen integrals in Ref. [83]. The method has been applied to polylogarithms in Ref. [93]. In the physical literature it dates back to Debye's work on the specific heat [94] in 1912. It is important to separate the cuts of the polylogarithms, which are either located on the positive or negative real axis. In deriving representation also mixed terms appear. Moreover, inã Here b k,l (x) andb k,l (x) are low order polynomials in ln(x), ln(1 − x) and Li 2 (x) with rational coefficients in x. For the function Li 2 (x) one uses the well-known Bernoulli-representation [90,93]. The functions r (i) k,l,m (x) are rational in x and c, a k , b k , d k , e k , f k ,ã k ,b k andd k are constants. In part of the region the above representation yields an even higher accuracy than double precision. The polynomial representations (5.46-5.51) may be further compactified using Horner's method [95] and are well suited to efficiently generate grids for further numerical use. This also applies to the corresponding OME and Wilson coefficient. All expressions have been derived using the package HarmonicSums. Numerical checks were performed using the code of Ref. [96] and the code HPL 2.0 [97].
We now consider the limiting behaviour of a (3),PS Qq (x) for large and small values of the momentum fraction x. In the limit x → 1 one obtains in the large-x region is illustrated in Figure 3. One should note that the factors of the various expansion coefficients are partly rather different, which gives preference to less singular terms. x For small values x, the function a   The asymptotic behaviour at small values of x is depicted in Figure 4. The leading term ∝ ln(x)/x does nowhere describe a (3),PS Qq (x), which is a common experience in many small-x studies, cf. [98,99]. An important term in the small x region is the 'next dominant' one ∝ 1/x, which has firstly been calculated in this paper 7 . The first two terms give a sufficient description up to x ≃ 5 · 10 −4 . At larger values, up to x ≃ 2 · 10 −2 in the small-x region all corrections ∝ ln k (x), k = 5...1 are required.
Let us reconsider the term ∝ ln(x)/x in a (3),PS Qq in view of the leading order small-x approximation. We mention that a rigorous theory of QCD in the small x limit has still not been worked out. This concerns, in particular, non-leading terms. As has already been known from the most singular terms ∝ 1/x, contributing to the pole 1/(N − 1) in Mellin-N space, in the unpolarized leading order splitting functions P gg (x) [100,101] and P gq (x) [102,103], they are related by the ratio C A /C F . In the ladder-approach using a physical gauge, this is the effect of exchanging one gluon-cell by that of a massless quark-cell [104]. In this way, leading small-x contributions to the constant part of the unrenormalized OMEs A Qg and A PS Qq might be related. This is indeed the case, which is easily seen for the expressions at 2-loop order in Mellin space [16], see also [11] for the expressions in x-space. The leading small-x contribution to the massive Wilson coefficient H 2,g has been calculated in [105]. Using the representation of the asymptotic heavy flavor Wilson coefficient given in Ref. [12] and expanding the result of [105] in the limit Q 2 ≫ m 2 using the corresponding perturbative representations of the BFKL resummed anomalous dimension [106,107], following e.g. [98], one obtains a prediction for the term of O(ln(x)/x) in a (3) Qg . This has been performed in [108]. Multiplying this expression by C F /C A agrees with the small-x limit of the exact expression (5.41), which is a new result of the present calculation. As has been shown above, the knowledge of this term is not enough for a quantitative prediction even in the small-x region and various more terms have to be computed. Following a method, developed by T. van Ritbergen [109] 8 , the authors of [108] used the fixed Mellin moments at 3-loop order calculated by Bierenbaum, Klein and one of the present authors in Ref. [12] using some sets of functions to determine a band of possibilities for values of a  To what extend these ansätze are compelling is hard to say, since the functional structure in a higher order calculation is difficult to predict and only understood after the calculation has been performed. In particular, the simplification of the initially large amount of generalized harmonic 8 Other approximate methods based on orthogonal polynomials to reconstruct x-shapes from moments have been known and were widely used even earlier, see Refs. [110]. We would like to mention that the determination of all these quantities, which are recurrent in N can be determined exactly knowing a finite number of moments as has been shown in [79]. The corresponding number of moments is, however, for mathematical reasons far larger than N = 6 as used in [108]. polylogarithms to the functions H andH is far from obvious. This applies also to the general functional structure. Furthermore, the former estimate [108] has been based on the unproven hypotheses of C F /C A scaling and the assumption that the leading small-x term for the Wilson coefficient can be determined in the way being described above. In Figure 5 we compare the range suggested for a (3),PS Qq in [108] with the exact result. In the important region of small values of x the guess in [108] has an uncertainty of ∼ 20%. At large values of x the grey area in Figure 5 winds around the exact result. This is enforced by known Mellin moments of Ref. [12] used to solve a linear system for the particular function set selected.
We now turn to the OME A PS Qq in the on-shell scheme for the quark mass. It receives contributions from O(a 2 s ) onward. We denote the logarithmic scale dependence by The matrix element reads in Mellin-N space as follows (5.55) The polynomials P i are given by P 34 = 5N 4 + 22N 3 + 49N 2 + 32N + 4 (5.56)       The difference between the OME in the MS-scheme and the on-shell scheme in N -space is given by The corresponding expression in x-space reads Here we have set the heavy quark mass equal in both schemes to obtain a more compact expression. The corresponding representations of the quark masses are given in [111]. The analytic continuation of the expressions in Mellin-N space can be obtained using the asymptotic expressions, which can be derived in analytic form [32,80], and the recurrence relations. Alternatively, the Mellin-inversion can be performed analytically and one may work in x-space. The corresponding relations for the OME are given in Appendix C.
In the numerical representation of the following figures the harmonic polylogarithms were calculated using the code HPLOG5 [90]. In Figure 6 the massive OME A PS Qq is shown as a function of x and Q 2 in the range of lower values of x. Here and in the following we illustrate the corrections to O(a 2 s ) and up to O(a 3 s ) using next-to-next-to-leading order (NNLO) parton distribution functions and the value of a s at NNLO, i.e. the term O(a 2 s ) refers not to the next-to-leading order (NLO) correction, but to the O(a 2 s ) in the NNLO correction. We refer to the value of a s (µ 2 ) as given in the parameterization [4] using the corresponding LHAPDF library. While in the small-x region the O(a 2 s ) term is negative, the NNLO correction is positive. Here and also in case of the heavy flavor Wilson coefficient in Section 6, one has to apply a Mellin convolution with the singlet quark densities at NNLO for a prediction in the variable flavor number scheme and in calculating the contribution to the structure function F 2 (x, Q 2 ) in the fixed flavor number scheme. Since the singlet distribution is decreasing towards larger values of x, the medium-x behaviour of the OME is important for the physical effect. Therefore, we enlarge the behaviour of the OME in the region of medium and large values of x in Figure 7. The O(a 2 s ) term remains negative and the NNLO correction turns to negative values approaching zero in the limit x → 1 from below. At larger values of x the NNLO correction becomes smaller than the O(a 2 s ) term. Unlike in the non-singlet case, we cannot present yet the matching in the VFNS to 3-loop order, since also the OME A Qg contributes here.

The Pure Singlet Wilson Coefficient in the Asymptotic Region
The pure singlet Wilson coefficient in the asymptotic region is given by Eq. (2.6). It also depends on the scale ratio Q 2 /µ 2 via In Mellin-N space it is given by 14) +320N + 144 (6.21) −23424N − 4032 (6.38) +259200. (6.41) In Appendix C we present the corresponding expression in x-space. Again one may also represent the structure function in Mellin-N space, using the corresponding evolution operators, cf. e.g. [98] and Wilson coefficient. This requests the representation of the corresponding (generalized) harmonic sums for N ∈ C, which can be realized using their analytic asymptotic representation to high accuracy and the known shift relations [32,80]. Figure 8 shows the x-and Q 2 -dependence of the pure singlet Wilson coefficient up to 2-and 3-loop order at Q 2 = µ 2 in the region of smaller values of x. While at low scales Q 2 ≃ 20GeV 2 the O(a 2 s ) term is positive in this region, it turns negative at higher scales. On the other hand, the corrections up to O(a 3 s ) are positive. At larger values of x, similar to the behaviour of the OME, also the corrections up to NNLO turn negative and, at even larger values of x undershoot the O(a 2 s ) contributions turning to zero for x → 0, as shown in Figure 9. The latter behaviour is of importance for the Mellin convolution with the quark-singlet distribution function. This is displayed in the contribution of the heavy flavor pure singlet contribution to the structure function F 2 (x, Q 2 ) in the case of the single heavy quark contributions at O(a 2 s ) and up to O(a 3 s ) in Figures 10 and 11, respectively, again setting µ 2 = Q 2 and referring to the onshell masses 9 . We apply the same setting for a s as in the case of the OME A PS Qq . Both corrections are negative. At lower scales Q 2 , the O(a 2 s ) effects are larger than those at NNLO, which give larger negative corrections at higher scales. To quote a few numbers, we obtain using the parton distribution functions [4]  range at HERA, x 0 = Q 2 /(sy) for y = 1, s ≃ 10 5 GeV 2 , one obtains the following corrections : −0.013 (x 0 = 0.0001), −0.026 (x 0 = 0.005), −0.023 (x 0 = 0.01), −0.019 (x 0 = 0.02), −0.007 (x 0 = 0.05) and 0.003 (x 0 = 0.1) for the contribution due to charm. In the whole kinematic region at HERA up to Q 2 ≃ 10 4 GeV 2 the bottom quark contribution is about one order of magnitude smaller than that by the charm quark.

Conclusions
We have calculated the O(a 3 s ) heavy flavor contributions to the flavor pure singlet OME A Qq , contributing to the matching relations in the VFNS and the corresponding heavy flavor Wilson coefficient for single heavy quark flavors in the asymptotic region Q 2 ≫ m 2 . As a by-product of the calculation we computed the complete pure singlet anomalous dimension at 3-loop order in an independent way for the first time. Our result agrees with the result in Ref. [87]. On the technical side of the calculation, we used the integration by parts package Reduze2 to reduce the Feynman integrals carrying local operator insertions to master integrals. The master integrals were calculated using different techniques. Most notably, we used differential equations, turned them into difference equations and solved using the packages Sigma, EvaluateMultiSums, SumProduction, HarmonicSums and OreSys. Part of these packages were also used to compute the nested sums obtained using representations by generalized hypergeometric functions and using Mellin-Barnes techniques. The massive OME and Wilson coefficient depend on new functions, which did not occur in the related 3-loop results having been dealt with previously in Refs. [23,24]. In the present result generalized harmonic sums contribute to the final expression in N -space. In x-space their effect manifests in harmonic polylogarithms of argument 1 − 2x. We studied the behaviour of the constant part of the unrenormalized massive OME, a Qq . This quantity has been newly calculated beyond the terms contributing to A (3) Qq and H PS q,2 implied by renormalization [12]. The leading small-x contribution of O(a 3 s ln(x)/x) is not describing this quantity, not even at values of x in the LHC region. To get physical values in the region of x ≃ 10 −4 one has to add the term O(a 3 s /x), which we also obtained as a by-product of the present calculation. To describe the region of larger values of x ≃ 10 −2 quite a series of logarithmic corrections in x have to be known. The pure singlet corrections both to structure function F 2 (x, Q 2 ) are negative and are largest in the small x region. The corrections due to bottom quarks are about one magnitude smaller than those by the charm quarks. We presented all quantities both in N and x-space for the use in deep-inelastic data analyses. The relations presented in the present paper can be obtained in computer-readable form on request via e-mail to Johannes.Bluemlein@desy.de.

A Integral families
In this appendix, we show the integral families that were implemented in order to perform IBP reductions for the operator matrix elements in the pure singlet case.
We give the corresponding set of propagators and depict the different topologies that are covered by each integral family, although only a few of these topologies are actually related to A (3),PS Qq , the rest being related to other OMEs. In the same way as for the diagrams on the right-hand side of Figure 2, a given bilinear propagator is indicated by a large dot on the corresponding line, and since the direction of the momentum is important in these propagators, this is also indicated in the diagrams. A solid line in the diagrams represents a massive particle (heavy quark), while a dashed line is a massless particle (gluon or light quark).
The names we have given to the families are arbitrary, of course. There is, however, a rationale behind the names we have chosen. The family names shown here start with a B, which indicates that these families are associated with a Benz-like topology.
In addition, they are constructed such that they also cover ladder topologies. Future calculations will require also crossed-box (non-planar) topologies, which we have labeled with names starting with a C, and are not shown here. The number that comes after the B (or the C) labels different routings of the mass in the propagators, while the letter that comes after this number labels the choice of three artificial propagators.
Family B5c (A

B Integrals
In the following we present representations of the generalized harmonic sums occurring in the pure singlet OME as Mellin transforms and in intermediary steps of the calculation, partly with different support, of (generalized) harmonic polylogarithms :  On the expense of more complicated arguments they can be represented in terms of classical polylogarithms up to weight w = 3. This is generally expected for 3-letter alphabets, see [72].