The logarithmic contributions to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\alpha _s^3)$$\end{document}O(αs3) asymptotic massive Wilson coefficients and operator matrix elements in deeply inelastic scattering

We calculate the logarithmic contributions to the massive Wilson coefficients for deep-inelastic scattering in the asymptotic region \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2 \gg m^2$$\end{document}Q2≫m2 to 3-loop order in the fixed flavor number scheme and present the corresponding expressions for the massive operator matrix elements needed in the variable flavor number scheme. Explicit expressions are given in Mellin \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document}N-space.

In the kinematic region at HERA, where the twist-2 contributions to the deep-inelastic scattering (DIS) dominate, cf. [35], 3 , i.e., Q 2 /m 2 10, with m = m c the charmquark mass, it has been proven in Ref. [37] that the heavy flavor Wilson coefficients factorize into massive operator matrix elements (OMEs) and the massless Wilson coefficients. The massless Wilson coefficients for the structure function F 2 (x, Q 2 ) are known to 3-loop order [38]. In the region Q 2 m 2 , where Q 2 = −q 2 , with q the spacelike 4-momentum transfer and m the heavy quark mass, the power corrections O((m 2 /Q 2 ) k ), k ≥ 1 to the heavy quark structure functions become very small.
In Ref. [39], a series of fixed Mellin moments N up to N = 10, . . . , 14, depending on the respective transition, has been calculated for all the OMEs at 3-loop order. 4 Also the moments of the transition coefficients needed in the variable flavor scheme (VFNS) have been calculated. Here, the massive OMEs for given total spin N were mapped onto massive tadpoles which have been computed using MATAD [41].
In the present paper, we calculate the logarithmic contributions to a series of unpolarized massive Wilson coefficients in the asymptotic region Q 2 m 2 to 3-loop order and the massive OMEs needed in the VFNS. These include the logarithmic terms log(Q 2 /m 2 ) and those parts of the constant term, which are obtained through renormalization, cf. [39]. Since the results for the OMEs A (3) gq and the fla-vor non-singlet Wilson coefficients for the structure function F 2 (x, Q 2 ) have been calculated completely in the meantime and have been published in [42,43], the corresponding expressions are not presented here. Results given in [36] will be discussed in Sect. 4. In the following, we set the factorization and renormalization scales equal μ F = μ R ≡ μ and exhibit the log(m 2 /μ 2 ) dependence on the Wilson coefficients, besides their dependence on the virtuality Q 2 . The logarithmic contributions are determined by the lower order massive OMEs [44][45][46][47][48][49], the mass and coupling constant renormalization constants, and the anomalous dimensions [50][51][52], as has been worked out in Ref. [39]. For the structure function F L (x, Q 2 ), the asymptotic heavy flavor Wilson coefficients at O(α 3 s ) for tagged heavy flavor were calculated in [53]. Here, we present the result considering the inclusive final state for this structure function. In both cases, the asymptotic representation becomes effective only at much higher scales of Q 2 [37] if compared to those for F 2 (x, Q 2 ), however. We first choose the fixed flavor number scheme to express the heavy flavor contributions to the structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ). This scheme has to be considered as the genuine scheme in quantum field theoretic calculations since the initial states, the twist-2 massless partons, can, at least to a good approximation, be considered as LSZ-states. The representations in the VFNS can be obtained using the respective transition coefficients within the appropriate regions, where one single heavy quark flavor becomes effectively massless. Here, appropriate matching scales have to be applied, which vary depending on the observable considered, cf. [54].
The constant contributions of the un-renormalized OMEs have been calculated in two cases A PS qq,Q (N ) and A qg,Q (N ), in Ref. [55]. Here, we present both the renormalized OMEs and the corresponding massive Wilson coefficients, which contribute first at 2-and 3-loop order. We also derive numerical results for the Wilson coefficients. The quantities being presented in the present paper derive from OMEs which were computed in terms of generalized hypergeometric functions [56,57] and sums thereof, prior to the expansion in the dimensional variable ε = D−4, cf. [58][59][60][61][62]. Finally, they are represented in terms of nested sums over products of hypergeometric terms and harmonic sums, which can be calculated using modern summation techniques [63][64][65][66][67][68][69][70][71][72][73][74][75][76]. They are based on a refined difference field of [77] and generalize the summation paradigms presented in [78] to multi-summation. The results of this computation can be expressed in terms of nested harmonic sums [79,80]. The results in Mellin N -space can be continued to complex values of N as has been described in Refs. [58][59][60][81][82][83][84].
It is the aim of the present paper to provide a detailed documentation of formulae in N -space for all logarithmic contributions and the terms implied by renormalization contributing to the constant part of the heavy flavor Wilson coefficients of the structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ) and the massive OMEs needed in the variable flavor number scheme up to O(α 3 s ) and to compare to results in the literature. Here, we refer to a minimal representation, i.e., we use all the algebraic relations between the harmonic sums and the harmonic polylogarithms, respectively, leading to a minimal number of basic functions. Based on the known Mellin moments [39] we also perform numerical comparisons between the different contributions to the Wilson coefficients and massive OMEs at O(α 3 s ) referring to the parton distributions [21]. The paper is organized as follows. In Sect. 2, we briefly summarize the basic formalism. The Wilson coefficients L PS q, 2 and L S g,2 are discussed in Sect. 3. As they are known in complete form we also present numerical results. In Sect. 4, the logarithmic contributions to the Wilson coefficients H S g,2 and H PS q,2 are derived and discussed. The corresponding inclusive Wilson coefficients for the longitudinal structure function F L (x, Q 2 ) in the asymptotic region are presented in Sect. 5. In Sect. 6, we compare the different loop contributions to the massive Wilson coefficients and OMEs for a series of Mellin moments in terms of the dependence on the virtuality Q 2 . Section 7 contains the conclusions. In Appendix A, the massive OMEs dealt with in the present paper which are needed in the VFNS are given in Mellin N -space.

The heavy flavor Wilson coefficients in the asymptotic region
We consider the heavy flavor contributions to the inclusive unpolarized structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ) in deep-inelastic scattering, cf. [85,86], in case of single electro-weak gauge boson exchange at large virtualities Q 2 . At higher orders in the strong coupling constant these corrections receive both contributions from massive and massless partons in the hadronic final state, which is summed over completely. In the latter case, the heavy flavor corrections are also due to virtual contributions. We consider the situation in which the contributions to the twist-2 operators dominate in the Bjorken limit. Here, transverse momentum effects of the initial state do not contribute. In the present paper, we consider only heavy flavor contributions due to N F massless and one massive flavor of mass m. 5 The Wilson coefficients are calculable perturbatively and are denoted by C S,PS,NS i, (2,L) x, N F + 1, Here, x denotes the Bjorken variable, the index i refers to the respective initial state on-shell parton i = q, g being a quark or gluon, and S, PS, NS label the flavor singlet, pure-singlet and non-singlet contributions, respectively. The massless flavor contributions in (1) may be identified and separated in the Wilson coefficients into a purely light part C i, (2,L) , and a heavy part by: (2,L) x, N F + 1, (2,L) x, N F , Q 2 μ 2 +H S,PS i, (2,L) x, N F + 1, (2,L) x, N F + 1, The heavy flavor Wilson coefficients are defined by L i, j and H i, j , depending on whether the exchanged electro-weak gauge boson couples to a light (L) or heavy (H ) quark line. From this it follows that the light flavor Wilson coefficients C i, j depend on N F light flavors only, whereas H i, j and L i, j may contain light flavors in addition to the heavy quark, indicated by the argument N F + 1. L S q, (2,L) can be split into the flavor non-singlet and pure-singlet contributions L S q, (2,L) = L NS q, (2,L) and at O(a 2 s ) only the non-singlet term contributes. The puresinglet term emerges at 3-loop order.
The heavy quark contribution to the structure functions F (2,L) (x, Q 2 ) for one heavy quark of mass m and N F light flavors is then given by, cf. [44], in case of pure photon exchange 6 1 x F Q Q (2,L) (x, N F + 1, Q 2 , m 2 ) = N F k=1 e 2 k L NS q, (2,L) x, N F + 1, (2,L) x, N F + 1, (2,L) x, N F + 1, (2,L) x, N F + 1, (2,L) x, N F + 1, The meaning of the argument (N F + 1) in Eq. (4) in the massive Wilson coefficients shall be interpreted as N F massless 6 For the heavy flavor corrections in case of W ± -boson exchange up to O(α 2 s ) see [89][90][91][92]. and one massive flavor. The symbol ⊗ denotes the Mellin convolution, 7 The charges of the light quarks are denoted by e k and that of the heavy quark by e Q . The scale μ 2 is the factorization scale, and f k , f k , and G are the quark, anti-quark, flavor singlet and gluon distribution functions, with The representations in Mellin N -space are obtained performing the transformation An important part of the kinematic region in case of heavy flavor production in DIS is located at larger values of Q 2 , cf., e.g., [93,94]. As has been shown in Ref. [37], the heavy flavor Wilson coefficients H i, j , L i, j factorize in the limit Q 2 m 2 into massive operator matrix elements A ki and the massless Wilson coefficients C i, j , if one heavy quark flavor and N F light flavors are considered. The massive OMEs are processindependent quantities and contain all the mass dependence except for the power corrections ∝ (m 2 /Q 2 ) k , k ≥ 1. The process dependence is implied by the massless Wilson coefficients. This allows the analytic calculation of the NLO heavy flavor Wilson coefficients, [37,46]. Comparing these asymptotic expressions with the exact LO and NLO results obtained from Refs. [4][5][6][7][8] and [1][2][3], respectively, one finds that this approximation becomes valid in case of F Q Q 2 for Q 2 /m 2 10. These scales are sufficiently low and match with the region analyzed in deeply inelastic scattering for precision measurements. In case of F Q Q L , this approximation is only valid for Q 2 /m 2 800, [37]. The different behavior if compared to F 2 (x, Q 2 ) is due to the emergence of terms ∝ (m 2 /Q 2 ) ln(m 2 /Q 2 ), which only vanish slowly in the limit Q 2 /m 2 → ∞. Here, we present the inclusive corrections for the structure function itself. The renormalized massive OMEs depend on the ratio m 2 /μ 2 , while the scale ratio in the massless Wilson coefficients is μ 2 /Q 2 . The latter are pure functions of the momentum fraction z, or the Mellin variable N , if one sets μ 2 = Q 2 . The mass dependence on the heavy flavor Wilson coefficients in the asymptotic region 7 Note that the heavy flavor threshold in the limit Q 2 m 2 is again x and not x(1 + 4m 2 /Q 2 ), which is the case retaining also power corrections. derives from the un-renormalized massive OMEŝ applying mass, coupling constant, and operator renormalization, as well as mass factorization, cf. Ref. [39]. The renormalized massive OMEs obey then the general structure The subsequent calculations will be performed in the MS scheme for the coupling constant and the on-shell scheme for the heavy quark mass m. The transition to the scheme in which m is renormalized in the MS-scheme is described in Ref. [39]. The strong coupling constant is obtained as the perturbative solution of the equation to 3-loop order, where β k are the expansion coefficients of the massless QCD β-function and μ 2 denotes the renormalization scale. For simplicity, we identify the factorization (μ F ) and renormalization (μ R ) scales from now on. In the subsequent sections, we present explicit expressions of the asymptotic heavy flavor Wilson coefficients in Mellin Nspace. They depend on the logarithms where μ ≡ μ F = μ R . Besides the Wilson coefficients (2), the massive OMEs are important themselves to establish the matching conditions in the variable flavor number scheme in describing the process of a single massive quark becoming massless 8 at large enough scales μ 2 , [39,44]. Here, the PDFs for N F + 1 massless quarks are related to the former N F massless quarks process independently. The corresponding relations to 3-loop order read, cf. also [44] 9 : The meaning of the symbolsf andf applied to a function here and in the following was defined in [39], Eqs. (2.3) and (2.4). Here, the N F -dependence of the OMEs is understood as functional and μ 2 denotes the matching scale, which for the heavy-to-light transitions is normally much larger than the mass scale m 2 , [54]. We will present the corresponding OMEs in Appendix A. The results of the calculations given in the subsequent sections have been obtained making mutual use of the packages HarmonicSums [95][96][97][98] and Sigma [63][64][65][66][67][68][69][70].
3 The Wilson coefficients L PS q,2 and L S

g,2
The constant parts of the un-renormalized massive OMEs A PS qq,Q and A S qg,A have been calculated in Ref. [55]. The OMEs contribute for the first time at 3-and 2-loop order, respectively, and stem from processes in which the virtual electro-weak gauge boson couples to a massless quark. Here and in the following N F denotes the number of massless flavors.
As a shorthand notation, we define the functioñ denoting the kinetic part of the leading order anomalous dimension γ (0) qg separating off the corresponding color factor. In Mellin N -space the Wilson coefficient L PS q,2 reads: with the polynomials For the massless 3-loop Wilson coefficients C k i, j we refer to Ref. [38]. Here and in the following, their expression will be kept symbolically. The harmonic sums are defined recursively by, cf. [79,80], The color factors are Likewise the Wilson coefficient L S g,2 is given by: where In all the representations of the massive Wilson coefficients and OMEs in N -space we apply algebraic reduction [99]. The 2-loop term in (28) is purely multiplicative and induced by renormalization only, while the 3-loop contributions require the calculation of massive OMEs. The above Wilson coefficients depend on the harmonic sums apart from those defining the massless 3-loop Wilson coefficients [38]. 10 In the above ζ l = ∞ k=1 1/k l , l ∈ N, l ≥ 2 denote the Riemann ζ -values, which are convergent harmonic sums in the limit N → ∞. In the constant part of the other Wilson coefficients it is expected that more complicated multiple zeta values emerge, which have been dealt with in [101].
In Eq. (28), denominator terms ∝ 1/(N − 2) occur. They cancel in the complete expression and the rightmost singularity is located at N = 1 as expected for this Wilson coefficient. Let us now consider both the small-and large-x dominant terms for both Wilson coefficients. Those of the massless parts were given in [38] before. Both Wilson coefficients contain terms ∝ 1/(N − 1). For simplicity, we consider the choice of scale Q 2 = μ 2 here. The expansion of the heavy flavor contribution, subtracting the massless 3-loop Wilson coefficients, denoted by L i , around N = 1(x → 0) and in the limit N → ∞(x → 1), setting Q 2 = μ 2 , is given by The corresponding limits for the contributions of the massless Wilson coefficients behave like 10 For the algebraically reduced representations see [100]. Q 2 =1000 GeV 2 Q 2 =100 GeV 2 Q 2 =20 GeV 2 Fig. 1 The O(a 3 s ) contribution by L PS q,2 to the structure function F 2 (x, Q 2 ) or m c = 1.59 GeV and the parton distribution functions [21] whereN = N exp(γ E ) and γ E denotes the Euler-Mascheroni constant.
In Fig. 1, we illustrate the contribution of L PS q,2 to the structure function F 2 (x, Q 2 ) using the PDFs of Ref. [21], cf. Eq. (4). Likewise, Figs. 2 and 3 show the corresponding contributions by L S g,2 at O(a 2 s ) and O(a 3 s ), respectively. The numerical results have been obtained using the representations in z-space, cf. Sect. 7, using numerical representations given in [102,103]. 11 Note that the O(a 2 s )-terms, cf. also Ref. [49], are smaller than those at O(a 3 s ), which is caused by terms ∝ 1/z in the momentum fraction z-space representation of the 3-loop contribution to L S g,2 , which are absent at 2-loop order.
These contributions emerging at 2-and 3-loop level are minor compared to the values of the structure function F 2 (x, Q 2 ), for which typical values are given in Table 1. 12 A global comparison of all heavy flavor contributions up to 3loop order can presently only be performed using the known number of Mellin moments, cf. [39], given in Sect. 6.
While the expression for L PS q,2 is the same in the MS and onmass-shell scheme to O(a 3 s ), L S g,2 , in its 3-loop contribution, 11 For other implementations see [104]. 12 Note that the kinematic region at small x probed at HERA is limited to values x ≥ Q 2 /(yS), with S 10 5 GeV 2 and y ∈ [0, 1].  Table 1 Values of the structure function F 2 (x, Q 2 ) in the low x region using the PDF-parameterization [21] x Q 2 = 20 GeV 2 Q 2 = 100 GeV 2 Q 2 = 1000 GeV 2 setting Q 2 = μ 2 . Here, we have identified the logarithms L M in both schemes symbolically. In applications, either the on-shell or the MS mass has to be used here. The Wilson coefficient H S g,2 , except for the constant contribution a (3),0 Qg , is given by: with the polynomials Our result for H S g,2 differs from the one given in Eq. (B.7) in z-space in Ref. [36], which misses the term given here in N -space. Note that the result of [36] is based on the calculation carried out by part of the present authors in Ref. [39], including the renormalization formulae being derived there. We have checked, however, that our result Eq. (72) is in full agreement with Eq. (2.15) of Ref. [39] and the moments having been calculated. For the logarithmic contributions to the pure-singlet Wilson coefficient H PS q,2 , leaving the constant part of the unrenormalized OME a PS,(3),0 Qq symbolic, we agree with the result given in [36]. Therefore, we refrain from presenting the corresponding expressions, also for the OME, contributing to the VFNS. The results on A (3),PS Qq will be given elsewhere.

The asymptotic Wilson coefficients for the longitudinal structure function
The Wilson coefficients in the asymptotic region have been calculated in Ref. [53] for exclusive heavy flavor production, consisting of three contributions only. In total also here five Wilson coefficients contribute and the expressions are modified in the inclusive case of the complete structure function F L (x, Q 2 ) cf. [39]. In Mellin N -space they read: with with and with As has been outlined for the 2-loop results in Ref. [37] already, the scales at which the asymptotic expressions are dominating are estimated to be Q 2 /m 2 800. They are far outside the kinematic region in which the structure function F L (x, Q 2 ) can presently be measured in deep-inelastic scattering. At LHeC [105], the region of large values of Q 2 can be probed and the corresponding expressions may be used in phenomenological analyses there.

Comparison of Mellin moments for the Wilson coefficients and OMEs
To compare the relative impact of the different Wilson coefficients on the structure function F 2 (x, Q 2 ), we will consider the Mellin moments for N = 2 to 10 in the following, convolved with the moments of the respective parton distribution functions in the flavor singlet case, i.e., the gluon G(x, Q 2 ) and quark-singlet density (x, Q 2 ) for N F = 3 and characteristic values of Q 2 . Since only a series of Mellin moments has been calculated at large momentum transfer Q 2 in Ref. [39], a detailed numerical comparison is only possible in this way at the moment. The numerical results for the moments of the contributing parton densities are given in Table 2. Note that for N ≥ 2, the moments for the singlet distribution are mostly larger than those of the gluon. We apply these parton densities to study the relative contributions of the different Wilson coefficients, normalizing to H S g,2 within the respective order in a s using the following ratios: In the numerical examples, we set e Q = e c = 2/3 and use m c = 1.59 GeV. Before we discuss quantitative results, a remark on the contributions by the color factor d abc d abc to the massless Wilson coefficients contributing to Eq. (2) is in order. Here we refer to Refs. [106][107][108] and [38,109]. For SU (n), one obtains It emerges weighted by 1/N c and 1/N A for external quark and gluon lines, respectively, with N c = n and N A = n 2 − 1.
In Refs. [106][107][108], this group-theoretic expression has been used, while in [38,109], a factor of 16 has been taken out and was absorbed into the Lorentz structure of the corresponding contribution to the Wilson coefficient. We agree with the N F -dependence as given in Refs. [106][107][108]. Furthermore, we note a typographical error in Eq. (4.13) of [38]. Here, the corresponding term reads correctly 13 Also in the pure-singlet case, the massless Wilson coefficients contain terms ∝ d abc d abc , although with a generally different charge-weight factor, cf. [106][107][108]. Let us now consider the relative impact of the individual massive Wilson coefficients. The ratios at O(a 2 s ) and O(a 3 s ) for different values of Q 2 and the moments N = 2 to 10 are given in Table 3. One first notes that at low values of Q 2 the ratio R(L S g,2 , H S g,2 ) changes sign, which is also the case for R(H PS q,2 , H S g,2 ) in the whole region up to Q 2 = 1000 GeV 2 . At O(a 2 s ) L S g,2 is small for low moments and grows to 24 % for N = 10 compared to H S g,2 at Q 2 = 100 GeV 2 , with lower In the case of the comparison of the massive OMEs, we normalize to A Qg with PDFs according to their appearance in the singlet and gluon transitions from N F → N F + 1 massless flavors in the variable flavor number scheme, cf. Eqs. (13)(14)(15): Here, we have included also the moments of the OME A NS qq,Q since it contributes to the singlet term in the VFNS. These ratios describe the relative impact, within the corresponding order in a s , of the massive OMEs in the variable flavor number scheme for the flavor singlet contributions. The numerical values for different scales of Q 2 are given in Tables 4 and 5 here. Its relative impact rises with Q 2 and amounts from ∼4 to 370 % for the R-ratio considering the lower moments N = 2 and N = 4 only.
Right after having obtained a series of moments for the massive OMEs at 3-loops in [39], it became clear that the logarithmic contributions are of comparable order to the constant term. Moreover, there is a strong functional dependence w.r.t. N , as displayed in Tables 2, 3, 4 and 5. To obtain a definite answer, the calculation of the constant parts of the un-renormalized OMEs a (3),0 i j as a function of N ∈ C is necessary. In particular, predictions for the range of small values x 10 −4 appear to be rather difficult using constraints by a low number of moments only.

Conclusions
We have derived the contributions of the massive Wilson coefficients to the structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ) in deep-inelastic scattering and the corresponding massive OMEs to 3-loop order in the asymptotic region Q 2 m 2 in Mellin N -space, except for the constant parts a (3),0 i j of the un-renormalized OMEs, which are not known for all quantities yet. Here, we retained both the scale dependence due to the virtuality Q 2 and the factorization and renormalization scales μ 2 , which were set equal. This allows for scale variation studies in applications. Two of the Wilson coefficients, L PS q,2 and L S g,2 , are known in complete form, and the corresponding results for L NS q,2 are given in [43]. In the variable flavor number scheme being applied to describe the process through which an initially massive quark transmutes into a massless one at high momentum scales, moreover, the matching coefficients A i j are needed. Here, A PS qq,Q and A qg,Q are known in complete form to 3-loop order and the results for A gq and A NS qq,Q are given in [42] and [43], respectively.
We have given numerical results for the Wilson coefficients L PS q,2 and L S g,2 . Using the available Mellin moments we have performed a numerical comparison of the different Wilson coefficients and operator matrix elements within the respective order in the coupling constant for the moments N = 2-10 and in the Q 2 range between 20 and 1000 GeV 2 . Here we have included also the OME from [42] and the Wilson coefficient and OME in the pure-singlet case for completeness. For numerical results in the non-singlet case see [43]. While some of the quantities studied are of minor importance, several others for the Wilson coefficients and OMEs are of similar size, which is varying in the kinematic range of experimental interest for present and future precision measurements. Even in case of the charm-quark contributions, the logarithmic terms are not dominant over the constant contributions in wide kinematic ranges, as, e.g., at HERA. Moreover, corrections of O (0.5-1 %) are important, given that the experimental resolution at HERA reaches 1 %. Much higher accuracies are expected at planned facilities as the EIC [33,34] and in part of the kinematical region at LHeC [105].
We have also derived the z-space expressions for all quantities given in the present paper, which were obtained from the N -space expressions using the package HarmonicSums [95][96][97][98]. Here, we have also reduced all expressions algebraically to obtain a compact basis representation [99]. These and the expressions in N -space may be obtained in form of computer-readable files on request via e-mail to Johannes.Bluemlein@desy.de.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited. Funded by SCOAP 3 / License Version CC BY 4.0.

Appendix A: The massive operator matrix elements in N-space
In this appendix, we present the massive OMEs in Mellin space to be used in the matching coefficients in the variable flavor number scheme Eqs. (12)(13)(14)(15). Thus far, the OMEs A PS qq,Q and A qg,Q are known completely. The other OMEs have been presented except for the 3-loop constant part a are presented elsewhere [42,43]. The transition matrix elements A PS qq,Q and A qg,Q are given by: with P 222 = 8N 7 + 37N 6 + 83N 5 + 85N 4 and with the polynomials Next, we present the OMEs, which are known except for the constant term in the un-renormalized massive OME at 3-loop order, a ,0 i j . The OME A Qg except for the term a (3),0 Qg reads: where P 320 = 333N 11 + 2331N 10 + 6556N 9 The operator matrix element A gg,Q except for the term a