The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_2(x,Q^2)$ and Transversity

We calculate the massive flavor non-singlet Wilson coefficient for the heavy flavor contributions to the structure function $F_2(x,Q^2)$ in the asymptotic region $Q^2 \gg m^2$ and the associated operator matrix element $A_{qq,Q}^{(3), \rm NS}(N)$ to 3-loop order in Quantum Chromodynamics at general values of the Mellin variable $N$. This matrix element is associated to the vector current and axial vector current for the even and the odd moments $N$, respectively. We also calculate the corresponding operator matrix elements for transversity, compute the contributions to the 3-loop anomalous dimensions to $O(N_F)$ and compare to results in the literature. The 3-loop matching of the flavor non-singlet distribution in the variable flavor number scheme is derived. All results can be expressed in terms of nested harmonic sums in $N$ space and harmonic polylogarithms in $x$-space. Numerical results are presented for the non-singlet charm quark contribution to $F_2(x,Q^2)$.


Introduction
The heavy flavor corrections to the structure functions in unpolarized deep-inelastic scattering yield large contributions in particular in the range of small values of the Bjorken variable x. Due to the current precision of the world deep-inelastic data which amounts for the structure function F 2 (x, Q 2 ) to O(1%) in a wide kinematic range, for the precision determination of the parton distributions [1], the strong coupling constant α s (M 2 Z ) [2] and of the mass of the charm quark m c [3], the heavy flavor corrections have to be known at 3-loop order. The heavy flavor corrections are known to 2-loop order in semi-analytic form [4]. 2 Due to a factorization of the heavy flavor Wilson coefficients into massive operator matrix elements (OMEs) and massless Wilson coefficients in the asymptotic region Q 2 m 2 [6], one may calculate the Wilson coefficients analytically in this region. In case of the structure function F 2 (x, Q 2 ) the asymptotic representation holds for scales Q 2 10 m 2 [6], while much larger scales are required for the structure function F L (x, Q 2 ). Analytic expressions for the heavy flavor Wilson coefficients in the asymptotic region and for the corresponding massive operator matrix elements were calculated in Refs. [6][7][8][9][10][11][12][13] at 2-loop order. The asymptotic 3-loop heavy flavor corrections to F L (x, Q 2 ) were computed in [14,15]. At 3-loop order, a series of Mellin moments has been calculated for the Wilson coefficients contributing to the structure function F 2 (x, Q 2 ) in [16] and for transversity in [17]. The logarithmic corrections and the parts of the constant contribution being determined by renormalization have been given in [15], cf. also [18]. For general values of the Mellin variable N , all OMEs corresponding to the color factors ∝ N F T 2 F C F,A were calculated in [19,20] and for the gluonic OME A gg,Q also for the terms ∝ T 2 F C F,A in [21]. Recently, also the OME A (3) gq (N ) has been computed [22].
In the present paper, we calculate the massive flavor non-singlet Wilson coefficient contributing to the structure function F 2 (x, Q 2 ) in the asymptotic region to 3-loop order and the associated massive operator matrix element A (3),NS qq,Q (N ). The latter quantity is given also for odd moments, which applies to the non-singlet contribution of the structure function g 1 (x, Q 2 ) as well as the charged current structure functions. We repeat the calculation for the tensor and pseudo-tensor operators which yield the corresponding OMEs for transversity [17,23]. As a by-product of the calculation, we obtain the contributions ∝ T F to the 3-loop anomalous dimensions in the vector and tensor case, in an independent calculation. In the case of transversity, this is the first calculation ab initio. We also present the heavy-to-light transition relations for the non-singlet distributions in the variable flavor number scheme (VFNS) to 3-loop order.
The paper is organized as follows. In Section 2, we give an outline on the flavor non-singlet OMEs and the massive Wilson coefficient to 3-loop order. In Section 3, we summarize technical aspects of the present calculation. The contributions to the 3-loop anomalous dimensions in the vector and transversity case are derived in Section 4. In Section 5, the OME in the vector case are given in Mellin space and in Section 6 the non-singlet Wilson coefficient is presented. Here we give also numerical results on its contribution to the structure function F 2 (x, Q 2 ). The nonsinglet transition functions in the VFNS are discussed in Section 7. The OME for transversity is calculated in Section 8. Section 9 contains the conclusions.
In Appendices A-C we present a series of master integrals used in the present calculation and summarize the corresponding representations in x-space.

Basic formalism
The massive flavor non-singlet operator matrix element is given by [16,17] A NS,MS qq,Q ( .p) N = q|O NS F,a;μ 1 ,...,μ n |q μ 1  Here denotes a light-like vector. The operators are contracted between massless quark states |q of momentum p. S denotes the symmetrization operator for the tensor indices, σ μν = (i/2)[γ μ γ ν − γ ν γ μ ], D μ is the covariant derivative, and λ r are the N 2 F − 1 matrices generating the SU (N F ) flavor group of N F massless quarks.
The massive flavor non-singlet Wilson coefficient in the asymptotic region Q 2 m 2 can be expressed in terms of the massive OMEs and the massless Wilson coefficients. The latter are known to 3-loop order [25]. The flavor non-singlet Wilson coefficient reads [6,16] L NS q,2 (N F + 1) = a 2 s A (2),NS qq,Q (N F + 1) +Ĉ (2),NS q,2 (N F ) + a 3 s A (3),NS qq,Q (N F + 1) + A (2),NS qq,Q (N F + 1)C (1), NS q,2 (N F + 1) +Ĉ (3),NS q,2 (N F ) . (2.7) Here the notation 'N F + 1' in the OMEs A qq,Q symbolically denotes that these are calculated at N F massless and one massive flavor. The 2-loop result has been calculated in Refs. [6,7] and the O(ε) term ā (2),NS qq,Q needed for the renormalization in [18]. The new contribution computed here is a (3),NS qq,Q both in the vector and transversity cases for the even moments. The odd moments refer to the corresponding polarized cases. The γ 5 -problem in the flavor non-singlet case is solved trivially since the operator is placed on the external line and the Ward-Takahashi identity [26] allows to map the vertex function into self-energy corrections, off the vertex.
The heavy flavor non-singlet contribution to the unpolarized deep-inelastic structure function F 2 (x, Q 2 ) for pure photon exchange is given by (2.8) where Here the diction 'flavor non-singlet term' for (2.8) was introduced in Ref. [10], rather meaning the individual flavor contribution. The non-singlet contribution in the sense of an associated operator would be f k +f k − Σ/N F , with Σ = N F l=1 (f l +f l ). For historical reasons we follow the former notion.
Before we present the results of the calculation, let us turn to its technical details first.

Diagrams and operator insertions
In order to calculate the operator matrix elements A (3),NS qq,Q and A (3),NS,TR qq,Q , we make use of the Feynman rules for operator insertions shown in Fig. 1, together with the standard Feynman rules of QCD, from which we construct and later evaluate the corresponding Feynman diagrams. There is a total of 112 diagrams needed for A (3),NS qq,Q , as well as for A (3),NS,TR qq,Q , which we generated using an extension of QGRAF [27] including also local operators [16]. The diagrams for A (3),NS qq,Q look exactly the same as the diagrams for A (3),NS,TR qq,Q , the only difference being in the explicit expressions for the operator insertions, where / is replaced by σ μν ν in the case of transversity. As we perform the renormalization for the reducible set of Feynman diagrams, cf. [16], the self-energy contribution has to be accounted for both the vector and transversity cases. Here ε = D − 4 is the parameter emerging in dimensional regularization. For the vector flavor non-singlet OME, it will lead to a vanishing 1st moment due to fermion-number conservation. A sample of the diagrams is shown in Fig. 2. The most complicated topologies are the Benzlike ones 4 of Figs. 2a-2c, with a massive triangle in the loop, and those with four massive lines as in the diagrams in Figs. 2d and 2e. No ladder or non-planar diagrams are involved in the calculation of these operator matrix elements, which simplifies matters considerably.
A FORM [29] program was written, cf. [16], in order to replace the propagators, vertices and operator insertions appearing in the output of QGRAF by the corresponding Feynman rules, and also to introduce the corresponding projectors and perform the Dirac matrix algebra in the numerator of the diagrams. After this, the diagrams end up being expressed as linear combinations of scalar integrals. For example, the diagram (also shown without labels in Fig. 2b) (3.2) is written as a linear combination of integrals of the form I a,b,c,N ν 1 ,...,ν 9 = Here m is the mass of the heavy quark and a.b denotes the Minkowski product. The denominators D 1 to D 8 are the inverse of the propagators appearing in the diagram. An additional auxiliary inverse propagator D 9 has been introduced in such a way that all possible scalar products of momenta k i .k j and k i .p (i, j = 1, 2, 3) can be uniquely expressed as linear combinations of all inverse propagators D 1 to D 9 . A scalar integral with irreducible numerators will be represented in this way by one or more negative integers among the exponents ν 1 to ν 9 in I a,b,c,N ν 1 ,...,ν 9 . In Eq. (3.3) there are also powers of the type .k i in the numerator. These come from the contraction of the / in the Feynman rule for the operator insertion in the diagram (3.2) with an internal momentum after the Dirac matrix algebra is performed in the numerator. At most one such power can appear in the case of this diagram, i.e., a, b, c ≥ 0 and a + b + c ≤ 1, but in the case of diagrams with 3-point and 4-point vertex operator insertions, higher powers can also appear. We will see in the next section that all operator insertions can be written in terms of artificial propagators, and introducing enough auxiliary propagators of this new type, the scalar products of internal momenta with will be expressed as linear combinations of artificial propagators.

Reduction to master integrals using integration by parts identities
We apply integration by parts identities [30] to the Feynman diagrams with operator insertions to reduce the problem to the calculation of a set of master integrals. The reductions were realized extending the C++-code Reduze 2 [31]. 5 In the following we describe the general implementation needed for the calculation of all the massive 3-loop operator matrix elements and focus later to the non-singlet case.

Mapping operator insertions to linear propagators
In order to systematically reduce loop integrals with operator insertions with Reduze 2, resummed contributions are introduced which effectively lead to linear propagators. Let us consider operator insertions involving fermion lines, where we suppress all factors independent of potential loop momenta. An insertion on a fermion line (FF) is mapped to scalar propagators introducing a new variable x and summing over N [34] according to 6 (3.5) The last equality defines a graphical representation for linear propagators. Similarly, we find for the second type of operator insertion (FFV) shown in Fig. 1 that the corresponding terms can be resummed into a generating function in the following way (3.6) after an appropriate shift in N and a partial fractioning. Note that the relative orientation of p 1 and p 2 in the linear propagators is important. We observe that it is possible to effectively decompose vertices with n edges into a combination of n − 2-vertices and n − 1 auxiliary edges. For n = 4, 5 we have (FFVV): and (FFVVV): Since we consider 3-loop diagrams, one of the four linear propagators in each term on the r.h.s. of the last mapping is redundant. In the corresponding diagrams, one of the legs will be external such that a combination of two linear propagators attached to the external leg can be split via partial fractioning.
For operator insertions at purely gluonic edges (VV) and vertices (VVV, VVVV, VVVVV) 7 we obtain similar patterns, although there are more combinatorial possibilities. In all cases, we find the relevant combinations of linear propagators represented by up to three connected edges in scalar diagrams with three-vertices only. The orientation of the linear momenta is such that 6 Here and in the following corresponding shifts in N have to be considered. 7 These terms will be required in future calculations of gluonic operators, see also Ref. [21]. they correspond to the direction of a transversal of these edges. If two consecutive edges differ only by an external momentum, one of them is ignored in the mapping to linear propagators.

Construction of integral families
In order to systematically handle our loop integrals, we index them using integral families. An integral family defines a set of scalar propagators which is complete and minimal in the sense that any scalar product of a loop momentum with a loop or external momentum can be expressed uniquely as a linear combination of inverse propagators, where the coefficients are kinematic invariants. As pointed out above already, this may require the introduction of auxiliary propagators, which are not fixed by the individual diagrams. In the following we describe how we construct a set of integral families such that all terms emerging from the diagrams can be indexed by powers of propagators from an integral family. Different routings of massive lines and the various combinations for the linear propagators require the introduction of many families. Note that in general one diagram will lead to a sum of terms which need to be matched to different integral families.
Our guiding principle is to cover all terms by a minimal number of families, which preferably have a large number of permutation symmetries. We start by considering diagrams without operator insertions. Possible 3-loop self-energy diagrams are described by Benz, ladder or crossed topologies (8 propagators) and their sub-topologies. They come in different variants, depending on the routing of the massive lines. Let us consider the following diagrams (3.9) The kinematics requires 9 propagators per family, that is, one auxiliary propagator. Requiring a large number of permutation symmetries and a small number of families, leads to the following choice of two families 8 -one which covers the depicted four planar topologies and one for the non-planar topology. We construct families for diagrams with insertions of type FFVVV, FFVV, FFV and FF for the massive fermion (in that order) and then for insertions of type VVVVV, VVVV, VVV, VV, where the latter scalar families also cover insertions involving massless fermions and ghosts.
Finally, we arrange the integral families such that families based on planar rather than nonplanar diagrams, with fewer massive propagators, or with higher permutation symmetries (in that order) are preferred. In this way we construct a total number of 24 integral families for the whole present 3-loop project (14 based on planar diagrams, 10 based on non-planar diagrams). It is curious to note that the above construction is quite stringent and leads to a relatively small set of integral families. As we will see, in the case of the OMEs A (3),NS,TR qq,Q only three families contribute.

Reduze 2: application and new features
We employ Reduze 2 [31] 9 for the reduction of loop integrals with integration-by-part identities and shift relations. Reduze 2 implements a distributed variant of Laporta's algorithm [38][39][40][41]. For this work, various aspects of the program were improved and extended. These improvements were implemented in a generic way and will, along with further new features, become publicly available with the upcoming Reduze release.
For this work, we extended various features of the program to support the linear propagators needed here. The previous version, Reduze 2.0, contained rudimentary support for bilinear propagators of the type 1/(q 1 .q 2 − m 2 ) where q 1 and q 2 are linear combinations of loop and external momenta. The new version of Reduze significantly extends support for this propagator type. Permutation symmetries are supported also if bilinear propagators are involved and explicitly verified to avoid user errors. In our example family B1a there is a symmetry under the following permutation of the propagators 1 ↔ 3, 2 ↔ 4, 6 ↔ 7, 10 ↔ 12, (3.11) corresponding to the shift of loop momenta k 1 ↔ k 2 . Restricting to permutations of subsets of propagators one can consider more general shifts. We extended the combinatorial shift finder (job setup_sector_mappings_alt) to determine such shifts of loop momenta also if bilinear denominators are permuted and/or a general crossing of external legs is involved. 10 These shifts are used to determine relations between sectors, possibly also between different families, and relations between integrals of specific sectors. In the present work, the only non-trivial crossing we need to consider is p → −p, which translates to x .p → −x .p at the level of the invariants. Furthermore, we improved the zero sector recognition for sectors with bilinear propagators. Applied to the families used in this work, these features reduce drastically the number of sectors for which a reduction needs to be performed.
In the presence of operator insertions, we find it more convenient to calculate integrals with higher denominator powers in comparison to integrals with additional numerators. We therefore choose an appropriate integral ordering to reflect this preference in the choice of the unreduced integrals in the reduction identities, that is, for the basis of master integrals. We added a new run-time option to the program which allows the user to select between various integral orderings. 9 For other public reduction programs see [35][36][37]. 10 Currently, we restrict the algorithms to permutations without additional minus factors, which is sufficient for the type of propagators considered in this work. Bilinear propagators with m 2 = 0, as present in different effective theories, can allow for mappings between propagators which involve additional minus signs: The new Reduze version provides a family finding algorithm, which systematically constructs families with a maximal number of permutation symmetries for a given set of propagators. Despite its current restriction to standard massless propagators, we found it useful for the construction of our families.
In the setup of the input files for Reduze, we implement the degenerate kinematics (3.10) via three incoming momenta p, −p, and appropriate rules for the scalar products. Some algorithms of Reduze assume non-degenerate kinematics. Certain modifications were required to allow Reduze to be used more flexibly for degenerate kinematics or other special setups. Most notably, this includes control over crossings of external legs to be used by the program via an explicit list.
Finally, the performance for long job queues and different usage aspects of the program were improved. We find it convenient to determine the master integrals needed in the calculation before the actual reduction is performed. This can be done more easily now using a new option we implemented in Reduze.
Let us consider the following example as an illustration for our IBP reduction procedure. In the case of the diagram (3.2), we can rewrite the integrals given in Eq. (3.3), which belong to family B1c, as 10 10 .
The operator insertion is resummed as If we appropriately introduce two additional linearly independent artificial denominators of this type, for example, then also the scalar products of with internal momenta in the numerator of Eq. (3.12) can be expressed as linear combinations of D 10 , D 11 and D 12 . In this way, all scalar integrals can be written as (3.15) A given scalar integral will be completely identified by the set of indices ν 1 to ν 12 . Of course, the specific form of the set of inverse propagators D 1 , . . . , D 12 will depend on the diagram we are considering. As outlined before, the auxiliary propagators in a given set should always be chosen so that the set is complete and minimal, which means that any scalar product of a loop momentum with , p or loop momenta can be expressed uniquely as a linear combination of the inverse propagators D 1 , . . . , D 12 defining a given integral family. We have found that all the diagrams needed for A (3),NS qq,Q and A (3),NS,TR qq,Q can be described by just three integral families, which are shown in Table 1.
In this representation, scalar integrals will be functions of x. A given integral I (x) will be written after reductions as a linear combination of master integrals J i (x) where the coefficients c i (x) are rational functions of x .p, the mass m and the dimension D.
Once an integral is obtained as a function of x, we can obtain the original integral we wanted to compute as a function of N by extracting the N th coefficient of the corresponding Taylor expansion in x, and doing the corresponding shift in N implied by Eqs. (3.5)-(3.8). It must be remarked that the coefficients c i (x) may contain poles in ε = D − 4, which will require the calculation of the corresponding master integrals beyond order ε 0 .
Since each diagram is written as a linear combination of scalar integrals I (x), the diagrams themselves have the structure of Eq. (3.16) in terms of master integrals, i.e., they are linear combinations of master integrals with rational coefficients in x. We can summarize the steps needed to calculate a diagram as a function of N as follows: Shift N as needed according to the corresponding operator insertion of the diagram.
In the next section, we will present some examples illustrating how the first step is performed. The remaining three steps are done by applying up-to-date computer algebra technologies. More precisely, after carrying out step 1, the master integrals I i (N ) are given in terms of definite multisums. Then in step 2 the obtained sums are simplified to expressions in terms of indefinite nested sums and products using the Mathematica package EvaluateMultiSums [42]. The backbone of this machinery relies on the package Sigma [43,44] that encodes advanced symbolic summation algorithms in the setting of difference fields [45][46][47][48][49][50][51][52][53]. In order to treat infinite summations these algorithms are supplemented by the package HarmonicSums [54,55] that can deal, e.g., with asymptotic expansions of the arising special functions.
After carrying out step 3, large expressions in terms of power series are given and a non-trivial task is to derive the N th coefficient in step 4. Here the package SumProduction [42] crunches the large amount of power series to a manageable number of basis sums and HarmonicSums calculates the N th coefficient of this compactified expression. The result is given in terms of definite sums which originate from the application of Cauchy-products. Therefore the packages Sigma and EvaluateMultiSums together with HarmonicSums are applied once more. As already mentioned earlier, we end up at expressions of the desired diagrams in terms of harmonic sums.
In Table 2, we show the list of all master integrals used for the reductions, where we identify each integral by giving the corresponding values of the indices ν 1 , . . . , ν 12 , see Eq. (3.15). We also indicate in the list the corresponding integral family and the order in ε to which these integrals need to be expanded. Of course, this list is not unique, and one is free to choose a different basis of master integrals. For convenience, we have chosen a basis of master integrals with no negative powers of inverse propagators.
Although the presence of the last three indices, ν 10 , ν 11 and ν 12 indicate in Table 2 that the integrals are functions of x, one may as well interpret them as functions of N by undoing the replacements as made in Eqs. (3.5)- (3.8). In Appendix A we give the results for the master integrals as functions of N .

Calculation of the master integrals
Compared with the master integrals required for the calculation of other operator matrix elements, the master integrals for A (3),NS qq,Q and A (3),NS,TR qq,Q listed in Table 2 are relatively simple. At most six propagators appear in these integrals (not counting the artificial propagators arising from operator insertions), and of these, at most three are massive, with the exception of the two integrals in family B5a, where four massive propagators appear. This simplifies the calculation of these integrals considerably and, in fact, many can be calculated solely in terms of Euler Betafunctions or a single sum of Beta-functions, which can then be performed using the package Sigma. Other more complicated cases required the solution of the master integrals in terms of hypergeometric functions [56][57][58], and the use of the corresponding series representations, in order to obtain a result in terms of a multiple sum to be solved by Sigma, EvaluateMulti-Sums and HarmonicSums. In this section we will show some examples describing the way these calculations are performed. To begin with, let us consider the integral J 2 . After introducing Feynman parameters and performing the momentum integrals, we obtain the following expression For simplicity, we have set the mass m and .p to 1, and an overall factor of i has been omitted. After doing a binomial expansion twice, integrating in w and in x in terms of a hypergeometric function, and using the following analytic continuation, Table 2 List of master integrals identified by the indices ν 1 to ν 12 according to Eq. (3.15), and by the integral families defined in Table 1. In the last column, we indicate the order in ε to which each integral needs to be expanded. An asterisk means that the integral is not only needed in the form presented in Eq. (3.15), but also in the form obtained by performing p → −p in this equation, which translates to x .p → −x .p at the level of the invariants.

Integral
Family B1b B1b At this point, we can use the series representation of the hypergeometric function since the argument is bounded between 0 and 1, and the series converges. We obtain . (3.20) This multiple sum is convergent and can be solved by our summation packages mentioned above in terms of harmonic sums. The result is given in Eq. (A.3).
Other integrals required the introduction of Appell hypergeometric functions. For example, after Feynman parameterization, the integral J 10 turned out to be given by where the ellipsis in the first line of the expression given above stand for the missing integrals in y and z 1 . Now we use the following analytic continuation, and obtain, It would be nice if we could use at this point the series expansion representation of the Appell hypergeometric function. However, if we do so here, we will get divergent sums, since the parameters of the F 1 function do not satisfy the conditions for convergence. We therefore rewrite our result using and get It is clear that the divergence of the F 1 function in Eq. (3.23) is partially due to the high negative power in the term (1 − z) −2− 3 2 ε appearing above. We can lower this power applying integration by parts in z, which leads to where These integrals will still produce divergent sums if they are re-expressed in terms of hypergeometric functions. In order to deal with K 2 , we now apply integration by parts in x, leading to where Now our original integral becomes, We see that integral K 1 does not need to be calculated, since it disappeared from the expression above. The integral K 2a can be done in terms of a 3 F 2 hypergeometric function, Integral K 3 is slightly more difficult. In order to do it, we first perform the following change of variables, which leaves the limits of integration from 0 to 1 unchanged, see also [59]. This leads to where we have dropped the primes in the variables y and z. The integral in x can be done in terms of a 2 F 1 hypergeometric function, which can then be expanded in terms of its series representation. Integrating then in the remaining variables, we obtain .
The sums in Eqs. (3.33) and (3.36) are convergent and can be performed using our summation toolbox, and after combining them according to Eq. (3.32) we obtain the result in Eq. (A.12). Except for integrals J 31 and J 34 , all of the other integrals were calculated in a similar way. The integral J 34 is, in fact, just a constant (independent of N ). In spite of that, it is still not easy to compute. 12 After Feynman parameterization we obtain, From here, we can obtain a Mellin-Barnes integral representation [61], using (3.38) in order to split the last term in Eq. (3.37). This allows to compute the Feynman parameter integrals in terms of Beta-functions. We get This integral can now be calculated with the help of the Mathematica package MB [62,63], which allows to analytically continue the integral for ε → 0 and afterwards perform the ε expansion. We obtain and 12 Integrals with fixed values of N or with no operator insertions can be obtained up to order ε 0 using the program Here N = N exp(γ E ) and γ E denotes the Euler-Mascheroni constant. The expression for b 2 is somewhat large, so it will be omitted here. These integrals can be computed by closing the contour to the left (or the right) and summing the residues. For example, in the case of b 0 , closing the contour to the left we get For more complicated expressions, the sum of residues can be performed using the package Sigma. The result for J 34 is given in Eq. (A.34). In the case of integral J 31 , we obtain We see that the N dependence factorizes outside of the integral, which is then just a constant. Except for this N -dependent factor, this integral is very similar to J 34 , and can be computed in an analogous way. The result is shown in Eq. (A.37).

The contributions to the 3-loop anomalous dimensions
In the calculation of the massive OME Eq. (2.5) one obtains from the contribution ∝ ln 2 (m 2 /μ 2 ) the 2-loop non-singlet anomalous dimension and from the term ∝ ln(m 2 /μ 2 ) the part of the 3-loop anomalous dimension ∝ N F . The anomalous dimensions γ ±NS have the representation where γ + NS (N ) is defined for even values of N and γ − NS (N ) for odd values of N . These are also the sets of values from which the analytic continuation to N ∈ C is performed. The anomalous dimensions are expressed in terms of harmonic sums [64] In the following, we drop the argument N of the harmonic sums and use the short-hand notation S a (N ) ≡ S a .

The vector case
We obtain the complete non-singlet anomalous dimensions at 2-loop order with the polynomials They agree with previous results given in Refs. [65]. For the contributions ∝ N F to the 3-loop non-singlet anomalous dimension the present calculation yields with the polynomials 14) The result agrees with the moments given in [66] and the general N result in [67].

Transversity
The complete 2-loop anomalous dimensions are obtained by .
They agree with the results obtained in Refs. [68]. For the contribution ∝ N F of the 3-loop anomalous dimensions we obtain γ (2) qq,NS, We agree with the moments given in [69,70] and note a typo in the 15th moment of [71]. In Ref. [71] the complete result for the 3-loop anomalous dimension was determined assuming the validity of maximal transcendentality studying the difference to the vector case, while the above result has been obtained in a direct calculation.

The flavor non-singlet massive operator matrix element: vector case
The new contribution beyond the terms determined by renormalization and factorization at 3-loop order is a (3),NS qq,Q , in Eq. (2.5). It is obtained as the constant part in the dimensional parameter ε of the unrenormalized 3-loop OME and given by Here and in the following expressions we will abbreviate some of the terms by the leading order non-singlet anomalous dimension normalized to the color factor C F . The constant B 4 reads which is related to a multiple zeta value, cf. [72], and the polynomials P i read 14) Here and in the following we reduce the contributing harmonic sums and harmonic polylogarithms algebraically [73]. Eq. (5.1) can be expressed in harmonic sums only reaching the level of weight w = 5. It is interesting to note that a part of the rational pre-factors rise ∝ N However, as can be seen form the asymptotic expansion 30) and the leading term behaves The flavor non-singlet OME in the on-shell scheme is given by   The corresponding expressions in x-space are given in Appendix B.
The difference between the OMEs in the MS-scheme and the on-shell scheme for the heavy quark mass in N -space is obtained by Here we used the same mass in both schemes symbolically to shorten the expression. The heavy quark masses in both schemes are given in Ref. [74]. In x-space it is given by Here H a (x) denote the harmonic polylogarithms [75] H b, a (x) =  H a (x). The +-prescription is defined by In the unpolarized case the leading term of the non-singlet OME in the small-x limit is obtained from Eq. (B.4) These terms correspond to the so-called leading poles at N = 0 in the flavor non-singlet case [76]. We finally derive the leading asymptotic behavior of the massive OME for large values of N , corresponding to the region x → 1. In the case of the on-shell scheme it is given by The corresponding expressions in the MS scheme are easily obtained using Eq. (5.53).

The massive flavor non-singlet Wilson coefficient in the asymptotic region
The heavy flavor unpolarized non-singlet Wilson coefficient in the region Q 2 m 2 , cf. Eq. (2.7), is given by The massless Wilson coefficient C (3),NS q,2 (N F ) has been given in [25] and the polynomials P i read     The Wilson coefficient depends on the logarithms where the renormalization and factorization scales μ F = μ R ≡ μ have been set equal. The corresponding expression in x-space is given in Appendix C in terms of harmonic polylogarithms. The analytic continuation can also be performed from the expressions in N -space directly, referring to their asymptotic representation and recurrence relations, cf. [77]. Let us comment on the pole-structure of Eq. (6.1). Because (6.1) is a flavor non-singlet quantity, one expects the rightmost pole to be situated at N = 0. Performing expansions around N = 2 and N = 1 one finds that the Wilson coefficient is finite there and the evanescent poles cancel.  To compare with results in the foregoing literature, we would like to consider the results of Ref. [6] to O(a 2 s ). Here the non-singlet heavy flavor Wilson coefficient was calculated in the case of tagged heavy quarks, given by the QCD-Compton process shown in Fig. 3. This scattering process is non-singular and can be calculated analytically for any ratio m 2 /Q 2 . The authors of [6] have given the following representation for Q 2 m 2 in terms of the massive OME and massless Wilson coefficients (1),NS qq (1),NS qq Note that in Eq. (6.1) the complete heavy flavor contributions were given, i.e. including also those due to virtual effects at 2-loop order. The difference term stems from the graphs shown in Fig. 4. In these graphs the final state contains massless partons only. However, the heavy flavor contributions to the deep-inelastic structure functions can only be defined consistently as the difference between the structure function involving all light and heavy quark contributions and an expression for the structure function considering only light flavors. The diagrams of Fig. 4 belong to the heavy flavor contributions at O(a 2 s ). In the asymptotic case μ 2 m 2 their contribution is given by −a 2 s β 0,Q ln 2,q . (6.37)  These terms exhibit a collinear singularity, being reflected in the μ-dependence of this contribution. Adding both terms yields Eq. (2.7) at 2-loop order.
In Fig. 5 we illustrate the contribution of the heavy flavor non-singlet Wilson coefficient Eq. (6.1) to the structure function F 2 (x, Q 2 ) accounting for the charm quark contributions at 2-and 3-loop order. The numerical calculation has been performed using the x-space representations given in Appendix C and the numerical representation of the harmonic polylogarithms given in [78]. 13 The contributions to O(a 2 s ) are significantly enhanced adding the O(a 3 s ) terms due to the gluonic and sea-quark contributions, despite the smallness of the additional factor a s . Overall, the NNLO corrections are below 1%. Note, however, that the present experimental precision for the structure function F 2 (x, Q 2 ) reaches O(1%). It will be even improved at high-luminosity facilities like the EIC [83] in the future.
In Fig. 6 we compare the predictions for the charm mass definition in the on-shell scheme and in the MS scheme. The corrections in the MS scheme lead to slightly larger absolute values than those in the on-shell scheme with very similar x-shapes.
The effect of retaining only the tagged heavy flavor contributions Eq. (6.36) is compared to the complete charm quark contributions to 2-loop order in Fig. 7. Here we set the factorization and renormalization scales to μ 2 = Q 2 . The absolute value of the correction is larger in the tagged case, compared to the inclusive corrections. We note that a separation of this kind at 3-loop order is in general not possible, without defining additional cuts, referring to the structure function approach as outlined in Ref. [16].

The variable flavor number scheme
The transition relation from N F → N F + 1 massless flavors of the flavor non-singlet distribution in the variable flavor number scheme is described by [16] 1) 13 For other implementations see [79][80][81]. where μ 2 denotes the matching scale. All the OMEs contributing to Eq. (7.1) are now known to 3-loop order, cf. also [15,19], which enables us to study its quantitative effects. The choice of the transition scale is process-dependent and usually significantly larger than the mass of the quark becoming light are requested to match the corresponding scattering cross sections, cf. [85]. This is obvious, because a heavy quark near production threshold is non-relativistic. In Fig. 8, we illustrate the flavor non-singlet, singlet and gluonic contributions to 3-loop order to the 4-flavor distribution x(u(x, μ 2 ) +ū(x, μ 2 )) as a function of x for μ 2 = 20 GeV 2 . The flavor non-singlet contributions are very small and show a weak x-dependence only, while the universal singlet and gluon distributions grow towards small values of x and are much larger. Due to this the distributions for the down and strange quarks are nearly the same.
In Fig. 9 the ratio of the momentum distribution function x(u +ū) is shown comparing the 4and 3-flavor scheme at NNLO for a variety of matching scales in the VFNS to O(a 2 s ) and O(a 3 s ). The ratio R(N F + 1, N F ) is defined by where the numerator is given by Eq. (7.1). At O(a 3 s ) the effects become larger, in particular at low matching scales and low values of x. In the region of larger x the predictions in both orders are closer. The overall effect is of O(±0.5%). Given the fact, that the structure function F 2 (x, Q 2 ) can be measured at the per-cent level, these effects are only somewhat below present experimental precision.
The corresponding ratio for the distribution x(d +d) is shown in Fig. 10. The results are very similar to those in Fig. 9, with a slightly larger effect at higher values of x. In comparison, the results for the distribution x(s +s) in Fig. 11

The flavor non-singlet massive OME: transversity
In a similar way to the non-singlet OME in the vector case, the OME for transversity is obtained referring to the operator Eq.
with the polynomials The transformation to the MS scheme is given by 20) in N -space. Again we have represented the mass in both schemes by the same symbol m. The corresponding expressions in x-space are given in Appendix B. Unfortunately the massless Wilson coefficients for transversity have not been calculated yet. Therefore the corresponding massive Wilson coefficients can only be given in the future, based on the present results.

Conclusions
The heavy flavor non-singlet corrections to 3-loop order in the asymptotic region Q 2 m 2 have been calculated in the vector case and for transversity as inclusive corrections in case of a single heavy quark at the time. We have presented the massive Wilson coefficient L NS 2,q (x, Q 2 ) and the transition matrix elements A qq,Q,NS(TR) occurring in the variable flavor scheme. They can be expressed in terms of harmonic sums in N -space and harmonic polylogarithms in x-space. At the technical side of the calculation 3-loop topologies up to massive Benz-type graphs had to be calculated. We used the integration-by-parts reduction for Feynman diagrams which carry local operator insertions, encoded in the package Reduze 2. The master integrals have been calculated using different methods like generalized hypergeometric functions, Mellin-Barnes techniques and differential equations, mapped onto summation problems. The latter ones were solved applying modern summation technologies encoded in the packages Sigma, Evalu-ateMultiSums, HarmonicSums and SumProduction.
As a by-product of the present calculation we obtained the complete non-singlet anomalous dimensions at 2-loop order and their contributions ∝ N F at 3-loop order, both in the vector and transversity cases. In the latter, we corrected a typo in the 15th moment in [71] and otherwise confirm earlier results in the literature. In case of transversity, the present result has been obtained in a first calculation ab initio.
The O(a 3 s ) corrections to F 2 (x, Q 2 ) amount to ∼ O(0.5%) in the lower x range in the 3-flavor scheme, with a small variance defining the charm mass in the on-shell or MS scheme. The present experimental accuracy for the structure function F 2 (x, Q 2 ) is about 1% in the small x region. We also discussed, at 2-loop order, differences in considering tagged heavy flavor vs. the inclusive heavy flavor corrections to F 2 (x, Q 2 ). The effect is smaller in the inclusive case. Finally, we considered the matching ratios in the VFNS going from 3 to 4 flavors, which now are available at O(a 3 s ) accuracy. In the ratio R(N F + 1, N F ) Eq. (7.2) the 3-loop effects are much larger than those at 2 loops in the small x region. At larger values of x the variation is broad and bounded by about ∼ ±1%. This applies in particular to low matching scales as applied in many phenomenological analyses.
The large formulae given in this publication are available in electronic form on request via e-mail to Johannes.Bluemlein@desy.de.

Appendix A. Master integrals
In the following we give the results for the master integrals as functions of the Mellin variable N . The integrals were calculated starting from N = 0, in other words, the shifts N − 1 → N and N − 2 → N were performed in the case of a fermion line insertion and a 3-point insertion, respectively. In some cases, the values at N = 0 are given separately, since the general N result diverges for this value of N . Each integral has been expanded up to the order in ε required for the present calculation. For simplicity, we have set the mass m and .p to 1, and the overall factor iS 3 ε has been omitted, where S ε is the spherical factor which emerges once at each loop order. At the end of the calculation in the MS renormalization scheme S ε is set to one. The different master integrals read , (A.2) 6(N + 1) , (A.6) , (A.7) The last expression diverges at N = 0. The correct result at this value of N is (A.12) , (A.14) This result is divergent at N = 0. Evaluating the integral separately at this value of N , we obtain

Appendix B. The operator matrix elements in x-space
The analytic continuation of the OME both in the flavor non-singlet vector case and for transversity leads to different expressions considering the even and odd moments from 3-loop order onwards, since these sets of moments correspond to different processes due to current crossing, cf. [86]. To account for this, we define the following functions