The three-loop polarized singlet anomalous dimensions from off-shell operator matrix elements

Future high luminosity polarized deep-inelastic scattering experiments will improve both the knowledge of the spin sub-structure of the nucleons and contribute further to the precision determination of the strong coupling constant, as well as, reveal currently yet unknown higher twist contributions in the polarized sector. For all these tasks to be performed, it is necessary to know the QCD leading twist scaling violations of the measured structure functions. Here an important ingredient consists in the polarized singlet anomalous dimensions and splitting functions in QCD. We recalculate these quantities to three-loop order in the M-scheme by using the traditional method of space-like off-shell massless operator matrix elements, being a gauge-dependent framework. Here one obtains the anomalous dimensions without referring to gravitational currents, needed when calculating them using the forward Compton amplitude. We also calculate the non-singlet splitting function ∆Pqq2,s,NS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta {P}_{\mathrm{qq}}^{(2),\mathrm{s},\mathrm{NS}} $$\end{document} and compare the final results to the literature, also including predictions for the region of small values of Bjorken x.


Introduction
Polarized deep-inelastic scattering allows to reveal the spin and angular momentum structure of nucleons [1]. At sufficiently high scales of the virtuality of the exchanged electroweak current Q 2 1 GeV 2 and for not too small or too large values of Bjorken x = Q 2 /(2p.q) [2] the twist-2 contributions dominate over the higher twist contributions [3,4] and target mass corrections [5,6]. Because of their strong variation, also the detailed control of the QED radiative corrections is required [7][8][9][10].
The scaling violations of the polarized deep-inelastic structure functions at the level of leading twist τ = 2 are determined by the anomalous dimensions of the composite quark and gluon operators [11][12][13][14] through the evolution of the polarized parton densities and the massless and massive Wilson coefficients. A major goal in measuring the polarized deepinelastic scattering process to high precision, e.g. at the EIC in the future [15] consist in the detailed measurement of the polarized parton distribution functions at a reference scale Q 2 0 together with the correlated measurement of the strong coupling constant a s (M 2 Z ) = α s (M 2 Z )/(4π) [16][17][18]. One of the important ingredients to this is the knowledge of the polarized singlet anomalous dimensions to highest possible order. The polarized nonsinglet anomalous dimensions are known to 3-loop order [19,20] and the polarized singlet anomalous dimensions were calculated to first [21][22][23][24][25], second [26][27][28][29] and third order in [30,31], using different computation techniques. In ref. [30] the calculation has been performed using the forward Compton amplitude, which required reference to gravitational subsidiary currents for the gluonic sector. A direct calculation of the contributions ∝ T F , 1 of the threeloop anomalous dimensions is possible using massive on-shell operator matrix elements [31], without reference to gravitational currents. However, this method only allows to calculate the anomalous dimensions ∆γ (2),PS qq and ∆γ (2) qg in complete form.
In the present paper we use the traditional approach of massless off-shell operator matrix elements (OMEs), allowing the direct calculation, as before for the non-singlet case JHEP01(2022)193 and the three-loop anomalous dimensions of transversity in ref. [20]. However, for this method the computational framework is gauge dependent. Still the anomalous dimensions in minimal subtraction schemes remain gauge invariant. Specific physical projectors have to be used, cf. also [29]. In the polarized case, unlike the unpolarized case [32], there is no mixing with the so-called alien operators [29,[32][33][34][35][36][37][38][39][40][41][42][43][44]. The main goal of the present paper is the re-calculation of all three-loop polarized singlet anomalous dimensions by a different method, being the first re-calculation of the polarized anomalous dimensions ∆γ (2) gq and ∆γ (2) gg . To some extent different steps in the computation technology are the same as in the non-singlet case, cf. [20] which will be referred to this paper and we only will describe new additional elements in the present paper. 2 At the phenomenological side we mention that, unlike the singlet case, the flavor non-singlet description of the structure functions g NS 1,2 (x, Q 2 ) is well explored to three-loop order, also including the heavy flavor contributions [45][46][47].
The paper is organized as follows. In section 2 we derive the structure of the physical part of the flavor singlet polarized unrenormalized off-shell OMEs to three-loop order. From their pole terms of O(1/ε) one can extract the singlet anomalous dimensions. We work in the Larin scheme [29,48] and perform finally the transformation to the M-scheme [29,30]. We also present the calculation details. Here a central method being applied is the method of arbitrary high moments [49]. In section 3 we calculate also the polarized non-singlet anomalous dimension ∆γ (2),s,NS qq , cf. also [50], and present the polarized singlet anomalous dimensions in section 4. Comparisons to the literature are given in section 5, including a discussion of the small x limit. Section 6 contains the conclusions and a new Feynman rule is presented in appendix A.

The unrenormalized polarized operator matrix elements and details of the calculation
The massless off-shell singlet OMEs are defined as expectation values of the local operators between quark (antiquark) ψ (ψ) and gluonic states F a µ 1 α of space-like momentum p, p 2 < 0. The OMEs are given by For the other definitions, here and in the following, we refer the reader to ref. [20].

JHEP01(2022)193
The Feynman rules of both QCD and of the local composite operators are given in [31,51,52] and have to be extended by that of the polarized local five-gluon vertex, given in appendix A.
We use the Larin scheme [48] 3 to describe γ 5 in D = 4 + ε-dimensions and express γ 5 by The Levi-Civita symbols are now contracted in D dimensions, The Larin scheme is one of the consistent calculation schemes in the polarized case. We will later transform to another scheme, the M-scheme. It is needless to say that the calculation of any observable requires to calculate also the Wilson coefficients in the same scheme, which implies that the extracted twist-2 parton density functions are obtained in this scheme, despite being universal w.r.t. the scattering processes in which they are used. Any analytic continuation of γ 5 or the Levi-Civita symbol to D dimensions violates Ward identities. One way to restore these consists in the calculation of scheme-invariant quantities. Let us add a few remarks on the M-scheme, cf. [29]. It is one of the consistent schemes for dealing with γ 5 in D dimensions. Historically, it has been motivated because here the supersymmetric relations hold to two-loop order [26][27][28][29], which are, however, not valid anymore at three-loop order, cf. (4.95) below and ref. [30]. The presentation of our results in the M-scheme has the reason to compare with the results of ref. [30]. Finally, one will perform the evolution of polarized structure functions in scheme-invariant form, as e.g. by studying the pair of observables g 1 (x, Q 2 ), ∂g 1 (x, Q 2 )/∂ ln(Q 2 ), cf. [69,70], or using one of the variety of the so-called DIS schemes, [71], cf. e.g. [72], to define the polarized parton densities, representing g 1 (x, Q 2 ) by purely quarkonic contributions, as at tree level.
The polarized operator matrix elements have the representation
In Mellin N space the Z-factor of a local singlet operator reads [52] (2.14) In (2.14) the terms ∆γ (2.15) The partly renormalized polarized singlet OMEs, ∆Ã phys ij , read The expansion coefficients a are given in ref. [88]. The anomalous dimensions in the M-scheme [29,30] with [29] z (1) and Let us now turn to the technical aspects of the present calculation, which has been widely automated and the corresponding chain of programs has been described to some extent in our previous paper [20]. A main point concerns the observation of the current

JHEP01(2022)193
crossing relations in the present case, cf. [11,89] and also ref. [88], projecting onto the contributions of the odd integer moments Through this one obtains a quadratic dependence on the resummation variable t. An expansion in t leads to the moments again. Our chain of programs, from the generation of Feynman diagrams to the final result, the polarized anomalous dimensions at three-loop order includes the packages QGRAF, FORM, Color, EvaluateMultiSums, Crusher, SolveCoupledSystems, Guess, Sage, Sigma, HarmoncisSums [52,. More details were given in ref. [20]. We determine the contributions to the anomalous dimensions due to the different color and multiple zeta values [122] individually. In the unfolding process also higher order terms a (k) ij in the dimensional parameter ε are needed to be calculated, cf. ref. [88]. We use the integrationby-parts relations [98,[124][125][126][127][128] to reduce the problem to master integrals. The method of arbitrary high moments [49] allows to find the needed difference equations, which are then solved using algorithms in difference ring theory [129][130][131][132][133][134][135][136][137][138][139][140][141]. To calculate initial values for the solution we use relations given in [142][143][144][145][146][147].
We turn now to some statistical characteristics of the present calculation. In the polarized flavor singlet case 125 irreducible diagrams contribute for ∆A  is extracted. The total number of irreducible diagrams in the polarized singlet case at three-loop order is larger by a factor of ∼ 22 than in the two-loop order. The reducible diagrams are accounted for by wave-function renormalization [84,86,87], decorating the OMEs at lower order in the coupling constant [29,88].
To calculate the anomalous dimensions ∆γ (2) ij we generated 3000 odd moments. It turns out that the determination of the largest recurrence requires 462 moments for ∆γ (2),PS qq , 989 moments for ∆γ (2) qg , 1035 moments for ∆γ (2) gq , 1568 moments for ∆γ (2) gg . In this way the anomalous dimensions are obtained. The largest difference equation contributing has order o = 16 and degree d = 304. These numbers are of the order obtained in the non-singlet case in ref. [101], where the largest difference equation contributing had order o = 16 and degree d = 192, and required 1079 moments. A moment based test-run for the polarized anomalous dimensions, like in the unpolarized case, has not been performed.
The overall computation time using the automated chain of codes described amounted to about 18 days of CPU time on Intel(R) Xeon(R) CPU E5-2643 v4 processors. Again we have only retained the first power of the gauge parameter to have a first check on the renormalization.
Because of the different requirements on the respective integrals, one distinguishes three contributions to the individual splitting functions in z-space which are defined in eqs. (38), (39) of ref. [20]. Also here we reduce the expressions to the algebraic basis, cf. [121], which has the advantage that only the minimal set has to be calculated in numerical applications [148][149][150]. We will not present the splitting function in explicit form, since the expressions are rather lengthy. They are given in computer-readable form in the supplementary material. The polarized singlet anomalous dimensions are given in section 4.

The polarized non-singlet anomalous dimension ∆γ (2),s,NS qq
In our previous paper [20] we had not yet calculated the non-singlet anomalous dimension ∆γ (2),s,NS qq , which emerges from three-loop order onward. It is best calculated using the vector-axialvector interference term in the forward Compton amplitude, corresponding to the associated structure function g − 5 (x, Q 2 ), ref. [89], corresponding to the difference in case of W − and W + charged current scattering. Due to its crossing relations it has even moments. This anomalous dimension has been calculated previously in ref. [50].
The corresponding gauge boson vertex is parameterized by and we consider the current interference term ∝ a · v. The projectors for the massless external quark lines of momentum p and boson lines corresponding to a tensor of rank two read The forward Compton amplitude depends on the invariants Q 2 = −q 2 and p.q = Q 2 /(2z) ≡ (Q 2 /2)y. The diagrams can be represented as formal power series in y cf. e.g. [20], section 2. The IBP-reduction in this case leads to three families and 101 master integrals in total. We consider the differential equations for the individual master integrals M k (y, ε), which are computed using the method described in [151]. After insertion of the master integrals into the amplitude and subsequent expansion in the dimensional parameter ε, the anomalous dimension can be determined by using the command GetMoment The corresponding splitting function reads Here ζ k denotes the Riemann ζ-function at integer argument k ≥ 2 and the color factor is normalized in the present case to d abc d abc /N c = 5/18 in QCD. The leading small z contribution of ∆P (2),s NS is ∝ ln 2 (z). This behaviour is, however, not dominant in kinematic regions being accessible at present. The asymptotic behaviour reaches the complete function up to 10% below z ∼ 10 −11 only. One more logarithmic order allows a description below z ∼ 10 −3 .

The polarized singlet anomalous dimensions
In the following we use the minimal representations in terms of the contributing harmonic sums and harmonic polylogarithms by applying the algebraic relations [121] between the harmonic sums and the harmonic polylogarithms. The polarized singlet anomalous dimensions can be represented by the following 23 harmonic sums up to weight w = 5 to three-loop order (4.1)

JHEP01(2022)193
This number is further reduced using also the structural relations [123,152] to at least 15 sums. In the present case the 10 sums suffice.
The splitting functions in z-space depend on the 26 harmonic polylogarithms 3) The harmonic sums are defined at the odd integers in the first place and the analytic continuation to N ∈ C is performed from there, [123,[153][154][155][156].
We obtain the following expressions for the polarized singlet anomalous dimensions in Mellin-N space, using the shorthand notation S a (N ) ≡ S a from one-to three-loop order.
Here we dropped the prefactor 1
The splitting functions in z-space are obtained by a Mellin inversion, cf. (2.35), and are given in computer readable form in the supplementary material.

Comparison to the literature
We confirm the results for the singlet anomalous dimensions calculated in [30], where the on-shell forward Compton amplitude has been used for the computation. The contributions ∝ T F have already been calculated independently as a by-product of the massive on-shell operator matrix elements in ref. [31], to which we also agree. The comparison to the large N F expansions of refs. [162,163] has been given in our previous paper [31] already, where all these terms are covered.
Let us finally consider the small z limit 5 of the present results and compare to the predictions given in [165,166]. In Mellin-N space these terms are given by the most singular contribution expanding around N = 0. One obtains at one-and two-loop order which agree with the evolution kernels given in refs. [165,166].

JHEP01(2022)193
At three-loop order we have in the M-scheme 3) and find a deviation both in case of ∆γ (2) qg and ∆γ (2) gq , already noticed in [167]. The relative deviation amounts to ∼ ±3.2%, with the difference of the expression in the M-scheme and the result obtained for the infrared evolution equation (IEE), δ∆γ (2),N →0 ij = ∆γ (2),N →0,M ij − ∆γ (2),N →0,IEE ij , δ∆γ (2),N →0 qg = 48 One may consider a different theory but QCD by setting C A = C F . In this case the prediction [165,166] for the evolution kernels agrees with the perturbative calculation of the anomalous dimensions. It has been the group of J. Kodaira [168], who also considered the effective Wilson coefficient in the non-singlet case ref. [169], 6 finding that these contributions are suppressed by a further power in N . This also applies to the effective Wilson coefficients in the nonsinglet case and therefore the evolution kernels of [169,170] agree. 7 By considering the expansion of the function to be interpreted as a matrix formulation for the effective Wilson coefficient in [165], The expansion coefficients F 0,k are 2 × 2 matrices only depending on color factors. Therefore, in the representation of [165], the Wilson coefficients are suppressed by one power in N , like in the non-singlet case, cf. [168,169]. To perform a full comparison one has to form two scheme invariant quantities. This is not really possible in pure polarized QCD, since there is only one structure function g 1 (x, Q 2 ) and also considering the Wilson coefficients [72,157]. 8 In [30] scheme invariant polarized evolution kernels for the contributions to the structure function g 1 and additional fictitious gravitational contributions in the gluonic channels have been formed for which the prediction in [165,166] hold. In data analyses or the phenomenological description of the polarized structure functions the consideration of only the leading small z terms in the evolution kernels is numerically not sufficient. Subleading terms dominate over the leading order effects, cf. [164,166,171]. 6 In using infrared evolution equations the authors of refs. [165,169,170] do not specify the factorization scheme, which complicates comparisons with calculations performed e.g. in the M-scheme. 7 After correcting ref. [170] in [171]. 8 At the level of twist-2 the structure function g2(x, Q 2 ) is not an independent quantity, but related by the Wandzura-Wilczek relation [172] to the structure function g1(x, Q 2 ). One might consider the physical pair {g1(x, Q 2 ), ∂g1(x, Q 2 )/∂ ln(Q 2 )}, cf. [69,70]. However, the so-called leading powers in N are here not of the same order.

JHEP01(2022)193 6 Conclusions
We have calculated the polarized three-loop singlet anomalous dimensions ∆γ (2),PS qq , ∆γ (2) qg , ∆γ (2) gq , ∆γ (2) gg , and the non-singlet anomalous dimension ∆γ (2),s,NS qq in Quantum Chromodynamics and agree with the results of refs. [30,31,50]. The singlet anomalous dimensions have been calculated by using the method of massless off-shell OMEs, which has been applied for this purpose for the first time. The calculation has been fully automated referring to the Larin scheme, performing a finite transformation to the M-scheme for the final results. Both schemes are valid to describe the scaling violations of the polarized structure functions, however, with a different outcome for the polarized parton distribution functions, which are scheme-dependent quantities. Comparing to predictions of the leading small z behaviour one finds deviations for the off-diagonal elements at three-loop order, from which one concludes that the calculation in ref. [165] is not in the M-scheme starting with three-loop order. To obtain the complete picture, one has to consider also the associated behaviour of the Wilson coefficients and to form scheme-invariant quantities. We also mention that partial checks on the polarized anomalous dimensions are possible from the pole terms of the single-and two-mass massive OMEs to three loop order, cf. [173][174][175][176][177][178][179][180]. Both the anomalous dimension and splitting functions are given in computer readable form in the supplementary material file ANOM3pol.m attached to this paper.
For the sums in the Feynman rule it is understood that the upper summation bound is larger or equal than the lower bound.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.