The polarized three-loop anomalous dimensions from on-shell massive operator matrix elements

We calculate all contributions ∝ T F to the polarized three–loop anomalous dimensions in the M–scheme using massive operator matrix elements and compare to results in the literature. This includes the complete anomalous dimensions γ ( 2 ), PS qq and γ ( 2 ) qg . We also obtain the complete two–loop polarized anomalous dimensions in an independent calculation. While for most of the anomalous dimensions the usual direct computation methods in Mellin N –space can be applied since all recurrences factorize at ﬁrst order, this is not the case for γ ( 2 ) qg . Due to the necessity of deeper expansions of the master integrals in the dimensional parameter ε = D − 4, we had to use the method of arbitrary high moments to eliminate elliptic contributions in intermediate steps. 4000 moments were generated to determine this anomalous dimension and 2640 moments turned out to be sufﬁcient.

and γ (2) qg . We also obtain the complete two-loop polarized anomalous dimensions in an independent calculation. While for most of the anomalous dimensions the usual direct computation methods in Mellin N-space can be applied since all recurrences factorize at first order, this is not the case for γ (2) qg . Due to the necessity of deeper expansions of the master integrals in the dimensional parameter ε = D − 4, we had to use the method of arbitrary high moments to eliminate elliptic contributions in intermediate steps. 4000 moments were generated to determine this anomalous dimension and 2640 moments turned out to be sufficient. As an aside, we also recalculate the contributions ∝ T F to the three-loop QCD β-function.

Introduction
The polarized three-loop anomalous dimensions γ (2) ij (N ) and splitting functions P (2) ij (z) govern the scale-evolution of the polarized parton distribution functions in Quantum Chromodynamics (QCD) at next-to-next-to-leading order (NNLO) and are of importance for precision predictions at epand hadron colliders, for the analysis of the different fixed target experiments, for the planned electron-ion collider EIC [1] and for RHIC. They are also instrumental for the measurement of the strong coupling constant α s (M Z ) [2] at these facilities and for the precise prediction of key processes like the polarized Drell-Yan process, jet production cross sections, and further processes. With the availability of the polarized 3-loop anomalous dimensions the present next-to-leading order data analyses of polarized deep-inelastic scattering data [3] can be promoted to the next-to-next-to-leading order level. Precision analyses of this kind are also relevant for the detailed study of the spin-composition of polarized nucleons (for reviews see [4,5]).
A first computation of the polarized three-loop splitting functions in the M-scheme was performed in Ref. [6]. The two-loop splitting functions have been known since 1995 [7,8]. In the flavor non-singlet case, the three-loop splitting functions P (2),NS− qq are the same as in the unpolarized case [9] and the contributions ∝ T F have been obtained in Ref. [10] as well. This also applies to transversity [10,11]. In the unpolarized case the three-loop splitting functions were calculated in Refs. [9,12] and all contributions ∝ T F were confirmed in independent massive calculations in Refs. [10,[13][14][15][16]. Already in 2010 we have computed the odd moments N = 1 − 7 of the polarized massive OME A (3) Qg and A (3) qg,Q , and more recently for N = 9. Before 2013 the corresponding moments for the massive OME A (3) gg,Q were calculated, as well as a similar number of the other OMEs at three-loop order. The whole set of moments remained unpublished, because an important detail in the definition of the massive OMEs with massless external quark lines first had to be understood.
In the present paper we compute the polarized three-loop splitting functions P (2),PS qq and P (2) qg and the parts ∝ T F of the three-loop splitting functions P (2) gq and P (2) gg from massive three-loop operator matrix elements (OMEs). They are necessary for the computation of the heavy flavor contributions to deep-inelastic scattering in the region of virtualities Q 2 much larger than the heavy quark mass squared m 2 . All splitting functions but P (2) qg are calculated by applying the techniques described in Refs. [17,18]. In the case of P (2) qg we use the method of arbitrarily high Mellin moments [19] to generate the moments of the O(1/ε)-pole of the corresponding OME, for which a recursion is obtained by using the guessing method [20]. This recurrence is finally solved by using the package Sigma [21,22]. The calculation is performed in the Larin-scheme [23]. As it turns out, in a massive calculation using on-shell massive operator matrix elements special care is necessary in treating massless external fermions, as we will explain later. At the end of the calculation we perform a finite renormalization to the M-scheme, cf. Ref. [24], to compare to the results in Ref. [6]. We would like to mention that the present calculation is thoroughly performed within QCD, while in Ref. [6] auxiliary graviton interactions had to be introduced to derive the gluonic anomalous dimensions γ (2) gq and γ (2) gg . The paper is organized as follows. In Section 2 we discuss the structure of the polarized unrenormalized three-loop massive OMEs, by which the three-loop anomalous dimensions can be calculated. If compared to the earlier literature, an important change has been necessary for the quarkonic projector, to obtain a consistent description of these quantities within the Larin scheme [23]. The calculation methods of the master integrals are summarized in Section 3. The finite renormalization from the Larin to the M-scheme [24] is described in Section 4. In Section 5 we present the polarized anomalous dimensions up to two loops, which can be obtained from the pole contributions of O(1/ε 3 ) and O(1/ε 2 ) of the massive OMEs. In Section 6 we present the contributions ∝ T F to the three-loop anomalous dimensions, which are the complete anomalous dimensions for γ (2),PS qq and γ (2) qg . In Section 7 we discuss the small z and large N F behavior of the splitting functions and anomalous dimensions and Section 8 contains the conclusions. In Appendix A we correct two of the operator Feynman rules given in Ref. [7]. Appendix B contains the splitting functions in z-space calculated in the present paper.

The polarized massive operator matrix elements
We will calculate the contributions ∝ T F to the three-loop anomalous dimensions using unrenormalized massive operator matrix elements. Their principal structure has been given in Mellin-N space in Ref. [25]. As an example we consider (2) gg,Q − N F a (2),PS Qq + a (2) Qg Qg . (1) Here we dropped the dependence on the Mellin variable N from all expressions for brevity. m denotes the bare heavy quark mass, ε = D − 4 the dimensional parameter, μ the factorization and renormalization scale, ζ l , l ∈ N, l ≥ 2 denote the values of the Riemann ζ function at integer argument, N F is the number of massless quark flavors, β i the expansion coefficients of the QCD β-function, β i,Q are related expansion coefficients associated to heavy quark effects, γ ij the expansion coefficients of the anomalous dimensions, and δm (l) k the expansion coefficients of the unrenormalized quark mass. The above quantities depend on the color factors C A = N C , C F = (N 2 C − 1)/(2N C ), T F = 1/2 for SU (N C ) and N C = 3 for QCD, cf. e.g. Ref. [25]. The coefficients a (k) ij denote the constant terms of the OMEs at k-loop order and ā (k) ij the corresponding terms at O(ε), cf. [26][27][28][29][30][31]. 1 In the constant terms O(ε 0 ) also multiple zeta values [34] contribute at fixed values of N . Furthermore, we use the convention The OME Â(3) Qg consists out of the one-particle irreducible part and a reducible part. These are related bŷ where ˆ (k) m 2 μ 2 denote the expansion coefficients of the massive contributions to the gluon vacuum polarization, see. [25]. The local operator insertions can be resummed into propagator-like structures, by virtue of an auxiliary parameter x, which form corresponding generating functions, cf. [10].
In total there are six massive OMEs, A qg,Q , A gq,Q and A (3) gg,Q , which we calculate in the polarized case. Furthermore, there is the non-singlet OME A (3),NS qq,Q which has already been calculated in Ref. [10]. The latter quantity, due to a known Ward identity, can be given in the MS scheme. From the poles O(1/ε 3 ) one obtains the one-loop anomalous dimensions and from the poles O(1/ε 2 ) the complete two-loop anomalous dimensions, while the contributions ∝ T F to the three-loop anomalous dimensions are extracted from the pole terms of O(1/ε).
To describe γ 5 in D = 4 + ε-dimensions, we work in the Larin scheme [23] 2 and express γ 5 by The Levi-Civita symbols are now contracted in D dimensions, In the calculation of the OMEs with external on-shell gluonic and quarkonic states we use the projectors P g and P q for the amplitudes Ĝ ab μν and Ĝ ij l . External ghost states do not contribute in the polarized states to three-loop order since the corresponding traces turn out to vanish. The gluonic projector is given by 1 The structure of the unrenormalized OMEs up to 3-loop orders given in Ref. [25] partly refers to earlier work by van Neerven et al. in the unpolarized case [32,33] to allow for a more direct comparison. This particularly concerns the definition of a (2) Qg and a (2) gg,Q . 2 For other schemes see Refs. [35]. For a discussion of the necessary finite renormalizations see [36].
Using it and (8) will imply the necessity of a finite renormalization of the two-loop anomalous dimensions γ (1),PS qq [37] and γ (1) gq [30], the structure and occurrence of which we did not understand for a series of years. Already in [27,37] and an early version of [29] we reanalyzed the pure singlet case at two-loop order to get agreement with the result of [26], which has been given there without presenting details. In the recent complete analytic calculation of the polarized massive two-loop Wilson coefficient in the whole kinematic range [28] it turned out that γ (1),PS qq did not receive a finite renormalization, as has been observed in the calculation of the polarized massless Wilson coefficient in [38][39][40] before. 3 For external quarkonic states a modified treatment compared to (10) therefore has to be applied. In the limit of a vanishing external light quark mass two bi-spinor structures survive, see [29]. They can be mapped to the following projector ε μνp tr p /γ μ γ νĜ ij l (11) which yields the proper definition in the Larin scheme in the case of the massive OMEs with massless external quark lines, unlike Eq. (10). The existence of a single projector (11) is of advantage since the calculation techniques described below only had to be modified minimally if compared to the unpolarized case. This also applies to the calculation of a series of fixed moments using MATAD [41]. The Feynman diagrams contributing to the massive OMEs were generated by the code QGRAF [42]. 4 The Dirac algebra has been performed using FORM [43] and the color configurations were calculated using the package Color [44]. The Feynman integrals were reduced to master integrals using the integration-by-parts (IBP) relations [45] implemented in the package Reduze 2 [46,47]. 5 There are different techniques available to calculate the master integrals, cf. Refs. [18,19], which we discuss in the next section.
Here a (2),NS qq,Q denotes the O(ε) contribution needed in the renormalization of the three-loop OME. We note that the use of the projectors (9, 11) for the massive OMEs allow to extract the contributions to the anomalous dimensions in the Larin-scheme 6 from all pole terms of O(ε −k ), k = 3, 2, 1. Their finite renormalization to translate to the M-scheme is described in Section 4.

The calculation methods
For the pole terms of the OMEs A (3),PS qq,Q , A (3),PS Qq , A (2) qg,Q , A (3) gq,Q and A (3) gg,Q the contributing master integrals can be calculated by the standard techniques such as the method of hypergeometric functions [50,51], the method of hyperlogarithms [52][53][54], the solution of ordinary differential equation systems [18,55,56] and the Almkvist-Zeilberger algorithm [57,58], being used in a combination, since in higher order in the dimensional parameter no elliptic integrals contribute. 7 Some of the simpler integrals have been calculated using Mellin-Barnes representations and using the codes in [59]. Most of the master integrals were already available from the calculation of the unpolarized three-loop anomalous dimensions in Ref. [16]. Only in a few cases some further differential equations had to be solved to obtain all master integrals. In all of the above methods corresponding sum representations have been derived which were solved using the difference-field techniques [60][61][62][63][64][65][66][67] of the packages Sigma [21,22], EvaluateMul-tiSums, SumProduction [68], and using HarmonicSums [58,[69][70][71][72][73][74]. This is, however, different for A (3) Qg . Due to the structure of the IBP relations some higher expansion in ε is necessary also to extract the term ∝ 1/ε. Here one would encounter elliptic contributions [75,76] by using the above techniques. We therefore apply the method of arbitrarily large moments [19] in this case. 8 Here one works in moment-space and the IBP relations are expressed in terms of recurrences for the master integrals. Using these relations one generates systematically higher and higher moments both for the master integrals and the operator matrix elements.
The projection to the analytic representation of the moments of the master integrals allows to treat also elliptic and higher structures. Finally one obtains the moments of the OME. They are used to derive a difference equation by the method of guessing [20] 9 implemented in Sage [79,80], based on very fast integer algorithms. We generated 2000 Mellin moments, which allowed to find most of the recurrences for all seventeen color-ζ projections. To determine the recurrences of the projections C F C A T F and C 2 A T F we used 4000 moments, out of which 2640 turned out to 6 The representation of the two-loop massive OMEs contributing to the structure function g 1 in the M-scheme are presented in Ref. [29]. 7 For a recent survey on these methods see [17]. 8 This method has been successfully applied also in a series of other calculations, cf. [16,77]. 9 For an early application to large systems in Quantum Field Theory see Ref. [78]. Table 1 Characteristics of the recurrences contributing to the anomalous dimension γ (2) qg . color/ζ order degree be sufficient. Here we refer to representations in terms of even and odd moments, with the even moments being unphysical. The analytic continuation is finally performed from the odd moments only. The characteristics of the recurrences for the different color-ζ factors contributing to the 1/ε term of the unrenormalized massive OME A (3) Qg are summarized in Table 1. For all the pole terms these recurrences are first-order factorizable and can be solved by applying the package Sigma. Here some color-ζ structures contribute for technical reasons, which cancel in the final expression.
All anomalous dimensions can be expressed by nested harmonic sums [69,70] To provide comparisons on a diagram-by-diagram basis we have calculated the first few Mellin moments for N = 1, 3, 5, 7, 9 using MATAD [41].

The finite renormalizations from the Larin to the M-scheme
We would like to compare to the results obtained in Ref. [6] which are given in the M-scheme. This scheme was defined in implicit form in Ref. [24]. Up to two-loop order it is the same as the one in which the results of Refs. [7,8] were obtained. The anomalous dimensions have the expansions in the non-singlet and singlet case At leading order, the anomalous dimensions are scheme-invariant. The finite renormalizations between the Larin and the M-scheme to three-loop order can be obtained following [24], see also [6], and are given by: with [24] z (1) and Specifically one obtains the following transformations: with the polynomials A priori, it has not been clear whether the use of the Larin scheme in the massive case leads to results which are equivalent to the HVBM scheme [35], which is known to occur in the massless case, cf. [6]. For the calculation of the anomalous dimensions it turns out that this is indeed the case. For the future it still remains to analyze all conditions implied by the Slavnov-Taylor identities, which are violated by dimensional regularization in both the Larin and HVBM schemes [23,35], in calculations of anomalous dimensions and Wilson coefficients from two-loop order onward.

The polarized anomalous dimensions up to two-loop order
Here and in the following we reduce the representations in Mellin-N space to bases by applying their algebraic relations, cf. [81]. We will also use the shorthand notation S a (N ) ≡ S a . The leading order anomalous dimensions are given by They are scheme-independent and agree with the results given in Refs. [82][83][84][85][86]. 10 The next-to-leading order anomalous dimensions are given by [6][7][8]28,29,31,39,[89][90][91][92] Here and in the following we only consider the non-singlet anomalous dimension γ The anomalous dimensions agree with results given in Refs. [6][7][8]28,39]. In the limit N → 1, which can be performed using HarmonicSums, we obtain with cf. [93][94][95][96]. Here, relation (62) is the consequence of the fact that the anomaly is maintained in its one-loop form, cf. [23].
As well the relation [6] γ (k),PS is verified. Furthermore, hold, where (66) is implied by fermion number conservation. The other first moments are given by

The contributions to the polarized three-loop anomalous dimensions ∝ T F
We will first present the anomalous dimensions in Mellin-N space for N an odd integer. Later this will form the basis to transform to the splitting functions to z-space. We consider the polarized massive OMEs A (2),PS qq,Q , A (2),PS Qq , A (2) qg,Q , A (2) Qg , A (2) gq,Q and A (2) gg,Q for odd values of N and derive the T F -dependent contributions to the anomalous dimensions from these quantities. After finally transforming from the Larin to the M-scheme one obtains: with andγ (2) with All anomalous dimensions agree with the results of Ref. [6].

The small z and large N F expansions
In the polarized case the so-called leading singularity of the anomalous dimensions is situated at N = 0 for the complete singlet matrix and in the non-singlet case, unlike the unpolarized case [99,100]. In Refs. [101][102][103] predictions were made for the small-z behavior of the flavor non-singlet splitting functions and in [103][104][105] and [106] for the polarized non-singlet and singlet splitting functions in QCD and QED, although not being fully clear within which scheme, using so-called infrared evolution equations.
Up to two loops, it turned out that the most singular contributions around N = 0 for the anomalous dimensions in the MS scheme are predicted. This is not the case at three loop order, where only the diagonal elements agree. 11 However, as shown in Ref. [6], if one interprets the predictions as physical evolution kernels, the 'leading' terms of the corresponding anomalous 11 A corresponding deviation has also been observed for the sub-leading BFKL anomalous dimensions in the unpolarized case [100]. dimensions agree. In the future it has still to be checked whether or not BFKL predictions will hold to even higher orders, cf. [107].
An interesting question concerns the numerical effect of the sub-leading corrections to the leading terms ∝ 1/N 5 , see Refs. [102,104,108]. One obtains The next sub-leading terms are more than canceling the leading terms, i.e. the leading terms alone yield basically nowhere dominant contributions in phenomenological applications numerically, despite the leading pole term of the complete calculation being correctly predicted. The finite renormalization between the Larin and the M-scheme does not affect the leading singular terms for the anomalous dimensions at N = 0.
We will now compare to predictions for the large N F terms given in Refs. [109,110]. Here we start with the flavor non-singlet case, [109], Eq. (3.5), as the generating function to understand the conventions used. The parameter μ is μ =D/2,D = 4 − 2ε,ε ≡ (4/3)T F ε and ε denotes the usual dimensional parameter. The function η o 1 is defined below in Eq. (3.12) [109]. The largest N F expansion coefficients to the non-singlet anomalous dimension are now obtained as the coefficient O(ε k+1 ), with the leading order term at O(ε). We agree with the corresponding non-singlet predictions.
The generating function for the coefficients is given in Eq. (5.3) [109], using the same definitions as in the non-singlet case above. Here γ (k) ij denotes the respective highest order term in N F . Expanding as in the non-singlet case, we reproduce all terms containing harmonic sums but not the rational term in Eq. (5.6) [109] after rescaling by a factor 64C F T 2 F as suggested by the expansion. Taking the rescaled result (5.6) [109] as it is and compare it to our three-loop result we obtain the difference accounting for the normalization of the anomalous dimensions in the present paper. This suggests that the finite renormalization to the M-scheme for the pure singlet term in Eq. (5.6) [109] still has to be done, while γ (2),NS qq and γ (2) gq , as suggested by Eq. (5.5) [109], are already in the M-scheme (note also a later discussion in Sect. 6 [109]). 12 We compare now the large N F term for the combination [110], to which we agree in the M-scheme. We finally would like to add a remark by J. Gracey: In respect of the large N F computation of [109], to accommodate the finite renormalization constant associated with the Larin scheme in perturbation theory, whereby chirality is recovered in strictly four dimensions, a correction had to be appended to the usual D-dimensional large N F critical exponent of the underlying operator. Such an additional piece is therefore by construction dependent on the procedure for handling γ 5 . Hence the prediction for perturbative coefficients beyond the leading one of [110] may not tally with those for an alternative γ 5 definition as appears to be the case here.

Conclusions
We have calculated the contributions ∝ T F to the polarized 3-loop anomalous dimension γ (2) ij (N ) and the associated splitting functions in a massive calculation, which is fully independent of the earlier computation in Ref. [6]. We agree with the previous results. To have the opportunity to deal with propagator-based representations only, we used formal Taylor series representations in terms of an auxiliary parameter x resumming the local operator insertions. As in the unpolarized case [16] before, we had to use the method of arbitrarily high moments [19] to deal with potential elliptic contributions in the necessary deeper expansions in the dimensional parameter ε in the case of the OME A (3) Qg . In the method of high Mellin moments [19] the moments are calculated recursively using the difference equation systems associated to the differential equations given by the IBP relations. Individual master integrals are only calculated in terms of moments. In all other contributions, standard techniques, cf. [17], are used in the calculation of the master integrals.
The universality of the QCD anomalous dimension allows to compute them within various setups. In the present calculation, they were obtained from the pole structure of massive polarized 12 Unlike for the terms at two-loop order, which were presented individually, only the combination (134) has been presented in [109]. Note that γ (1)PS qq is scheme independent, while γ (2)PS qq is not. OMEs. The calculation of these OMEs is part of an ongoing project with the final goal to compute the massive polarized Wilson coefficients for deep-inelastic scattering in the region Q 2 m 2 .
The anomalous dimensions and splitting functions presented in this paper are also given in Mathematica format in ancillary files to this paper. The conversion to maple or FORM-inputs [43] is straightforward.

Acknowledgement
We thank J. Gracey, F. Herren, S. Moch, D. Stöckinger, A. Vogt and S. Weinzierl for discussions. This work has been funded in part by EU TMR network SAGEX agreement No. 764850 (Marie Skłodowska-Curie) and COST action CA16201: Unraveling new physics at the LHC through the precision frontier. A. De Freitas has been funded in part by the Innovatives OÖ 2020 programme. The Feynman diagrams were drawn using Axodraw [112].

Appendix A. Feynman rules for local operator insertions
Two of the Feynman rules given in Ref. [7] had to be corrected. 13 They correspond to the vertices given in Fig. 1.
Here all momenta are ingoing. The quark-quark-gluon-gluon operator reads, with N ∈ N defining the Mellin moment, The generators of the color group are denoted by t c and g is the strong coupling constant with a s = g 2 /(4π) 2 13 We thank J. Smith for a corresponding communication related to Ref. [26] several years ago, with which we agree.
Here the symbols f abc denote the QCD structure constants. Nevertheless the two-loop anomalous dimensions calculated in [7] are correct.

Appendix B. The splitting functions
The splitting functions are related to the anomalous dimensions by the Mellin transform Since distribution-valued components are contributing, the Mellin transform is used as follows In the following we present the complete polarized splitting functions to O(α 2 s ) and the T F -dependent polarized splitting functions at O(α 3 s ). Here P (2),PS qq and P (2) qg are the complete splitting functions.
All splitting functions can be expressed by harmonic polylogarithms [111], which are given by the iterated integrals over the alphabet Again we reduce to the algebraic basis, cf. [81]. As a shorthand notation we use H a (z) ≡ H a .

B.1. The splitting functions to two-loop order
The leading order splitting functions are given by The next-to-leading order splitting functions read We obtain the following contributions ∝ T F to the next-to-next-to-leading order polarized splitting functions