The Three-Loop Splitting Functions in QCD: The Helicity-Dependent Case

We present the next-to-next-to-leading order (NNLO) contributions to the main splitting functions for the evolution of longitudinally polarized parton densities of hadrons in perturbative QCD. The quark-quark and gluon-quark splitting functions have been obtained by extending our previous all Mellin-N calculations to the structure function g_1 in electromagnetic deep-inelastic scattering (DIS). Their quark-gluon and gluon-gluon counterparts have been derived using third-order fixed-N calculations of structure functions in graviton-exchange DIS, relations to the unpolarized case and mathematical tools for systems of Diophantine equations. The NNLO corrections to the splitting functions are small outside the region of small momentum fractions x where they exhibit a large double-logarithmic enhancement, yet the corrections to the evolution of the parton densities can be unproblematic down to at least x about 10^{-4}.

A little more than ten years ago, we completed the calculation of the third-order (next-to-next-toleading order, NNLO) corrections P (2) ik , i, k = q, g, for the helicity-averaged (unpolarized) case [12,13].
These calculations were performed in the approach of Refs. [14,15] where physical quantities, specifically structure functions in inclusive deep-inelastic scattering (DIS), are calculated via forward amplitudes in dimensional regularization [16][17][18][19]. In order to access also the lower row of the NNLO flavour-singlet splitting-function matrix, i.e., P (2) gq and P (2) gg , in a third-order calculation, this procedure requires the inclusion of a process other than standard gauge-boson exchange DIS. The method of choice, cf. Ref. [20], was to include DIS via a scalar φ coupling directly only to gluons via φG μν a G a,μν , where G μν a is the gluon field strength tensor, as realized in the Standard Model by the Higgs boson in the limit of a heavy top quark and five massless flavours [21,22].
A corresponding calculation was performed six years ago for the structure function g 1 in polarized photon-exchange DIS, which is sufficient to extend the determination of the helicitydependent (polarized) splitting functions [23][24][25] to NNLO for the upper-row quantities P qq and P qg . Since we had no access to the corresponding lower-row splitting functions, these results were only briefly discussed in Ref. [26]. There is no helicity-sensitive analogue to the above Higgs-boson exchange in the Standard Model or an effective theory derived from it (initially a pseudoscalar χ with a χε μνρσ G μν a G ρσ a coupling to gluons was tried, which however cannot probe spin information either, as also χ is a scalar under the rotation group).
This leaves only working in supersymmetry, as in Ref. [27] for the determination of the NNLO quark-gluon antenna function, or considering DIS by graviton exchange. We have chosen to adopt the second option, which is easier to implement in our setup and offers additional information and checks by accessing all four splitting functions P ik , as well as their unpolarized counterparts, and a full set of physical evolution kernels for both the unpolarized and the polarized case.
The basic formalism for graviton-exchange DIS has been developed in Ref. [28]; for a recent application see also Ref. [29]. There are three structure functions H 1,2,3 in the unpolarized case, of which three combinations can be formed which are analogous to F 2 (no gluon contribution at order α 0 s ), F φ (no quark contribution at order α 0 s ) and F L (neither) in gauge-boson and scalar DIS. In the polarized case there are two structure functions, H 4 and H 6 , where H4 = H 4 − H 6 and H 6 involve only the quark and gluon distributions, respectively, at the leading order, in perfect analogy with the system (F 2 , F φ ) that we employed for obtaining the unpolarized splitting functions.
We have performed complete second-order calculations of all these quantities. At three loops, however, graviton exchange leads to a large number of integrals with a higher numerator complexity than encountered in the calculations for Refs. [12,13,26]. Hence repeating the step from fixed-N Mellin moments [14,15] to all-N results would require a lot of time and/or considerably improved algorithms. We have therefore resorted to calculating P (2) gq and P (2) gg for fixed (odd) values of N . Substantial improvement in our diagram handling and in the FORM [30][31][32] implementation of the MINCER program [33,34], see Ref. [35], together with the availability of sufficient computing resources, have enabled us to completely determine P (2) gq (N ) for 3 ≤ N ≤ 27 and P (2) gg (N ) for 3 ≤ N ≤ 25 (the N = 1 moments are not accessible in this calculation [28]), and both for specific colour factors up to N = 29.
Initially the extension to high moments was intended to facilitate approximate x-space results, analogous to but much more accurate than those obtained in Ref. [36] based on the moments of Ref. [37] for the unpolarized case, which would suffice at all x-values relevant to 'spin physics' in the foreseeable future. Similar to the somewhat simpler case of transversity in Ref. [38], however, it turned out that it is possible to reach values of N for which even the most complicated parts could be determined completely from the moments and additional endpoint information, in particular the suppression of P ik (x) − P ik (x) by two powers of (1 − x) in the threshold limit x → 1 in a suitable factorization scheme. The crucial step in this determination is the solution of systems of Diophantine equations for which we have, besides in-house tools coded in FORM, made use of a publicly available program [39] using the LLL-based [40] algorithm described in Ref. [41].
Consequently we are now in the position to present the complete NNLO contributions P (2) ik to the helicity-difference splitting functions in perturbative QCD. The remainder of this article is organized as follows: In Section 2 we set up our notations and discuss aspects of the second-order calculations and results relevant to our determination of the third-order corrections which we turn to in Section 3. Our N -space results for P (2) ik are presented in Section 4, and the corresponding x-space expressions in Section 5, where we also briefly illustrate the numerical size of the NNLO contributions to the evolution of polarized parton densities. We summarize our results in Section 6. Some additional information on scheme transformations and graviton-exchange DIS is collected in Appendices A-C. A brief account of this research has been presented before in Ref. [42].

Notations and second-order results
The unpolarized and polarized parton densities of a longitudinally polarized nucleon are given by and where f + i and f − i represent the number distributions of the parton type i with positive and negative helicity, respectively, in a nucleon with positive helicity. Here x denotes the fraction of the nucleon's momentum carried by the parton, and μ the mass-factorization scale which can be identified with the coupling-constant renormalization scale without loss of information.
The scale dependence of the quantities in Eqs. (2.1) and (2.2) is governed by the renormalization-group evolution equations d d ln μ 2 ( )f i x, μ 2 = ( )P ik α s μ 2 ⊗ ( )f k μ 2 (x) (2.3) where ⊗ stands for the Mellin convolution in the momentum variable, given by with a s ≡ α s (μ 2 ) 4π . (2.6) Using symmetries, the system (2. 3) of 2n f + 1 coupled integro-differential equations, where n f denotes the numbers of effectively massless flavours, can be reduced to 2n f − 1 scalar flavour non-singlet equations and the 2 × 2 system for the polarized gluon density f g (x, μ 2 ) and the flavour-singlet quark distribution The quark-quark splitting function P qq in Eq. (2.7) can be decomposed as P (n) qq (x) = P +(n) ns (x) + P (n) ps (x) (2.9) into non-singlet and pure singlet components. The former is related by P + ns = P − ns to an unpolarized quantity calculated in Ref. [12], the latter starts only at n = 1 and is specific to the present polarized case. It is often convenient to consider the Mellin transforms of all quantities, given by The complete next-to-leading order (NLO) contributions P (1) ik for the quantities in Eq. (2.7) have been derived almost 20 years ago in Ref. [23] in N -space using the OPE and in Refs. [24,25] in x-space, using the lightlike axial-gauge approach of Refs. [2,3]. Some years ago, we have checked these results, and obtained P (2) qq and P (2) qg , by extending the calculations for Refs. [13,43] to the structure function g 1 in polarized DIS which was first addressed beyond the first order in Ref. [44]. All these calculations used dimensional regularization, and thus needed to address the issue of the Dirac matrix γ 5 in D = 4 dimensions which enters via the quark helicity-difference projector.
The calculations in Ref. [23] used the 'reading-point' scheme for γ 5 [45]; those in Refs. [24,25] were carried out primarily with the 't Hooft/Veltman prescription [46,47], but included checks also using the so-called Larin scheme [48,49], / pγ 5,L = 1 6 ε μνρσ p μ γ ν γ ρ γ σ , (2.11) where the resulting contractions of two ε-tensors are evaluated in terms of the D-dimensional metric. All our calculations have been carried out using the Larin scheme which is equivalent to the 't Hooft/Veltman prescription for the present massless case. Quantities calculated using Eq. (2.11) need to be subjected to a factorization scheme transformation in order to arrive at expressions in the standard MS scheme [50,51], for example where we have switched to a matrix notation in N -space and suppressed all function arguments. Denoting the perturbative expansion of the transformation matrix by , (2.13) the transformation (2.12) of the coefficient functions C g 1 and the parton densities f leads to (2) , P (0) + Z (1) , P (1) L − Z (1) , P (0) Z (1) for the splitting functions in the MS scheme, where [a, b] denotes the standard matrix commutator. Here β 0 and β 1 are the leading two coefficients in the expansion of the beta function of QCD, s β , (2.15) which to NNLO is given by [4,5,[52][53][54][55] with C A = n c = 3 and C F = (n 2 c − 1)/(2n c ) = 4/3 in SU(n c = 3). β 0 and β 1 are schemeindependent in massless perturbative QCD; β 2 is given in the MS scheme adopted in this article.
The transformation matrix has been determined to NNLO in Ref. [56] as Its non-singlet entries can be fixed by the relation between the corresponding coefficient functions for g 1 and the structure function F 3 which is known to order α 3 s [57]; the critical part is the pure-singlet part for which, as far as we know, only that one calculation has been performed so far. For the convenience of the reader the results are included in Appendix A. For z (n) qg = z (n) gg = 0, Eq. (2.14) leads to the following transformations of the NLO and NNLO splitting functions: (1) gg,L + P (0) qg z (1) gq (2.18) and P (2) qq = P (2) qq,L + β 0 z (1) (2) gq , P (2) qg = P (2) qg,L + P (1) qg,L z (1) qq + P (0) qg z (2) qq , gg,L + P (1) qg,L z (1) gq + P (0) qg z (2) gq . (2.19) These expressions are reduced to the standard scheme transformation of Refs. [23][24][25]56] by dropping all contributions with z (1) gq or z (2) gq ; it will become clear below why these terms have been included in Eqs. (2.18) and (2.19).
It is instructive to consider the x → 1 threshold limit of the splitting functions. It is expected that the physical probability of a helicity flip is suppressed by two powers in (1 − x) in this limit [58]. Hence the differences (2.20) should be suppressed, in a 'physical' factorization scheme, by a factor of (1 − x) 2 , or 1/N 2 in N -space, relative to the respective sums which behave (modulo logarithms) as (1 − x) −1 or N 0 for ik = qq, gg and (1 − x) 0 or N −1 for ik = qg, gq. For the scheme-independent leading-order (LO) splitting functions, the differences (2.20) read The corresponding NLO results for the MS splitting functions [23][24][25] are given by Interestingly, as already noted in Ref. [26], all 10 terms in Eq. (2.23) can be removed by including the simple additional term z (1) gq in the NLO scheme transformation (2.18). The splitting functions P (1) qg (x) and P (1) gq (x) are shown, together with their unpolarized counterparts, in Fig. 1 in the standard scheme, from now on denoted by 'M' wherever required, that uses only Eq. (2.17) and an alternative scheme ('A') that also includes this additional term.
The issue of the physical large-x behaviour of the helicity-dependent quark-gluon splitting can be addressed by studying suitable flavour-singlet physical evolution kernels (or physical anomalous dimensions) for structure functions in unpolarized and polarized DIS. Gravitonexchange DIS, for which the basic formalism was worked out in Ref. [28], provides a sufficiently large set of structure functions. It is convenient to combine and normalize four of these functions as  [23][24][25] ('M') and after including an additional term z (1) gq in the transformation (2.14) from the Larin scheme ('A'), which removes all (1 − x) 0,1 terms from the quantity δ (1) gq (x) in Eq. (2.23).
with H2 = H 2 − 4H 3 in the unpolarized case, and with H4 = 2(H 4 − H 6 ) in the polarized case, where we have changed the x n prefactors relative to Eq. (31) of Ref. [28] such that (C u ) ij = (C p ) ij = δ ij at LO. The corresponding NLO coefficient functions can be found in Appendix B. The physical-kernel matrices K a , a = u, p (for the renormalization scale μ 2 R = Q 2 ) are obtained from the coefficient functions, the beta function (2.15) and the respective unpolarized (P u = P ) and polarized (P p = P ) splitting functions, cf. Eq. (2.7), by The expansion of this result to order a 3 s can be read off from Eq. (2.14) for Z = C a . We have performed complete two-loop calculations of these structure functions, recovering both the unpolarized and polarized NLO flavour-singlet splitting functions from gravitonexchange DIS, and used these results to obtain the NLO physical kernels K (1) u (x) and K (1) p (x). The respective off-diagonal elements for the systems (2.24) and (2.25) are compared in Fig. 2. It is clear, also from the corresponding analytical results, that also the large-x limits of the kernels K (1) 32 (x) and K (1) 64 (x) corresponding to the splitting functions ( )P (1) gq are consistent with the expectation of Ref. [58]; hence Eq. (2.23) is indeed a unphysical feature of the standard transformation to the MS scheme.

Determination of the third-order corrections
As before, we have calculated inclusive DIS via the optical theorem, which relates the probe (q)-parton (p) total cross sections (with Q 2 = −q 2 > 0 and p 2 = 0) to forward amplitudes, and a dispersion relation in x that provides the N -th moments from the coefficient of (2p · q) N [14,15]. For the splitting functions P (2) qq and P (2) qg we have extended the three-loop all-N calculations of Refs. [12,13] to the photon-exchange structure function g 1 . As discussed in Ref. [26], a large number of additional integrals, arising from a fairly small set of top-level integrals with higher numerator powers, had to be calculated for this extension; their determination took several months.
The situation is far worse in the case of graviton-exchange DIS, which is our means to access also P (2) gq and P (2) gg , in terms of both the complexity and the number of new top-level integrals. We have therefore not tried a direct all-N calculation in this case, but managed to set up a two-step procedure with the same result. The first step is a calculation of fixed-N moments for the structure functions in polarized graviton-exchange DIS, as in Refs. [14,15] using the MINCER program [33,34], but up to much higher moments in particular for H 6 , cf. Eq. (2.25). The second step is the determination of the all-N expressions for P (2) gq and P (2) gg from the moments calculated in the first step together with insight into the structure of these functions.
In order to drive the first step to a point where the second became possible, and its results could be verified by one or two yet higher moments, improvements had to be made in our diagram preparation and the MINCER code, see also Ref. [35]. The diagrams were generated, as before, with a special version of QGRAF [59]. Unlike in our previous calculations, however, the diagrams with the same group-invariant colour factor, the same topology and subtopology (see below), and the same flavour structure have been combined in the 'diagram' files which are managed, as before, using the database program MINOS [60]. In this way the number of thirdorder diagrams has been reduced from 5176 to 1142 and from 15 208 to 1249 for the quark and gluon contributions, respectively, to H 4 and H 6 . The combined diagrams take roughly as much time as the most difficult individual diagram in the set, which leads to an overall gain in speed by a factor of three to five.
The overall most demanding subtopology, in terms of execution time and required disk space, is NO 25 (see Fig. 3), i.e, the most difficult p-flow in the most difficult three-loop topology. Also notable are the LA 14 (also shown in Fig. 3), O4 57 , O2 26 cases, where the momentum p flows through four internal lines, and the three-line BE 57 and BE 28 'Benz' cases. The largest diagram calculated took about 10 7 CPU seconds and required 6.7 TB of disk space for the projection on N . The results for 3 ≤ N ≤ 25 were employed for obtaining the all-N expressions for P (2) gq and P (2) gg . For checking these expressions, the quark case was computed completely at N = 27 and in the 'planar limit' C A − 2C F → 0 at N = 29, and the gluon case for the C 3 A terms at N = 27 and N = 29. The latter was possible since most of the slowest and largest diagrams do not contribute to this colour factor, which is the most complicated one in terms of the structure of the splitting function.
Most of the diagram calculations were performed on the ulgqcd cluster in Liverpool, using TFORM [31,32] with 16 workers on more than 200 cores; the hardest diagrams at the highest values of N were calculated on a new high-end computer at NIKHEF. For the previous optimization of MINCER we were also able to use a multi-core workstation at DESY-Zeuthen.
As an example, we show the non-ζ 3 parts of the moments 3 ≤ N ≤ 25 of the C 3 F part of P In order to obtain, with certainty, the analytical forms of P (2) gq (N ) and P (2) gg (N ) from only 12 moments, we need to make use of additional constraints on the structure of these functions. At least up to NNLO, the splitting functions can be expressed in terms of harmonic sums [61], see also Ref. [62], which can be recursively defined by The sum of the absolute values of the indices m k defines the weight of the harmonic sum. Assigning a weight m to the un-summed denominators which can be expressed as differences of two harmonic sums of weight m, the N n LO splitting functions include terms up to weight 2n + 1. For example, the C 2 F n f contribution to P (2) in the standard MS scheme [56], where all harmonic sums are understood to be taken at argument N . Here we have also made used of the first of the abbreviations for the N -dependence of the lowest-order splitting functions, cf. Eq. (4.2) below.
If the unpolarized counterpart of Eq. (3.5) is written down in the same notation, the first two lines are the same except for the replacement of p qg by p qg = 2D 2 − 2D 1 + D 0 . The same holds for the C A C F n f and C 2 A n f contributions. As in other results in massless perturbative QCD, the number of harmonic sums is reduced by the absence of sums with index −1. This leaves Fig. 4. The NNLO differences δ (2) gq (N ) = P (2) gq (N ) − P (2) gq (N ) for the non-n f and n 1 f terms in the M and A schemes for C A = 3 and C F = 4/3, compared to the unpolarized result. The symbols show moments calculated using MINCER, the solid and dashed lines the exact all-N results presented below. As at NLO, cf. Fig. 1, the M-scheme difference turns negative at large N . seven sums of weight 3, of which one is missing in Eq. (3.5) but not the corresponding C A C F n f and C 2 A n f expressions. Half of their in principle 28 coefficients with D 0,1 and D 2 0,1 are fixed by the 1/N 2 suppression of the difference δ (2) qg in Eq. (2.20), which is found to hold separately for each harmonic sum. Taking into account the lower-weight sums, this large-N behaviour relates as many as 24 coefficients to the unpolarized result for each of the three non-n f colour factors.
Another crucial feature of Eq. (3.5) and all other available results for splitting functions is that all coefficients are integer in a suitable normalization. E.g., after eliminating all terms linear in D 0 and D 1 using the 1/N 3 large-N behaviour, the remaining coefficients in Eq. (3.5) are integers once factors of 2 w−3 have been bracketed out of the terms with sums of weight w < 3. Consequently the equations relating the remaining coefficients to fixed-N moments are Diophantine equations, and far less that n equations are required to determine n unknown coefficients. While there are a few additional constraints, on the coefficient of the D 5 0,1 and D 4 1 terms corresponding to the ln 5 x and x ln 5,4 x small-x logarithms and the remaining coefficients of S 1,1,1 (N ), see below, it is clear that it is vital for the determination of P (2) gg (N ) to have an extension of the A-scheme of Fig. 1 to NNLO, in order not to miss out on those 24 large-N constraints.
The double-logarithmic S 1,1,1 and S 1,1,1,1 contributions to P (2) gq (N ) can be derived from the calculations of polarized graviton-exchange DIS, without any reference to the unpolarized results, from the single-log threshold enhancement of the physical kernel K p in Eq. (2.26), cf. Ref. [63]. An additional scheme transformation that removes those contributions to δ gq is found to be The assumption that this remarkably simple transformation leads to is consistent with the results for N ≤ 25 is illustrated in Fig. 4 for the n 0 f and n 1 f contributions in QCD.
The physical kernels for the system (H4, H 6 ) also allow to settle another issue observed in Ref. [26], the apparent partial disagreement of the leading small-x logarithm of P (2) qg (x) with the old resummation result of Ref. [64]: the ln 4 x contribution to K4 6 agrees perfectly with that prediction, which clarifies its proper interpretation, see also Refs. [65][66][67]. Consequently it should be possible to use the prediction of Ref. [64], via K 64 , also for P (2) gq (x). Furthermore the x ln 5 x and x ln 4 x terms of this function can be fixed by extending the analysis of the small-x limits of the unfactorized expressions in Ref. [68] to the present case, see also Ref. [69].
Finally we need to briefly address the issue of denominators other than D 0 and D 1 , as occurring in the eighth line of Eq. (3.5), and with sums to weight 3 in its C A C F n f counterpart. Due to the different leading-order structure, there are far fewer such terms here than in the unpolarized case. Terms with D 2 in P (2) qg (N ), D −1 in P (2) gq (N ) and D −1 D 2 in P (2) gg (N ) do neither affect the prime-number decomposition of the denominators of the odd-N moments, e.g., the N = 17 moments do not involve a factor 1/19, cf. Eq. (3.1), nor can they lead to an overall pole at N = 1.
We are now ready to turn to the determination of the all-N expressions. The structure of the critical C 3 F , C A C 2 F and C 2 A C F parts of P (2) gq is analogous to Eq. (3.5) discussed in detail above. With the coefficients of the weight-4 sums fixed by the unpolarized result [13], we are left with 2 × 32 coefficients of sums at weight 3 and below combined with powers of D 0 and D 1 , recall Eq. (3.4), plus at most 11 sums combined with D −1 . The large-N suppression of δ (2) gq in the A-scheme and the other endpoint constraints fix 29 or 30 of these coefficients (depending whether or not D −1 S 1,1,1 (N ) is included in the basis set), leaving up to 45 unknown integer parameters.
We have developed FORM tools for analyzing the prime-number structure of the moments, see Eq. (3.1), and deriving relations between the remaining parameters using the Chinese remainder theorem [70]. These tools have proved sufficient, sometimes together with a brute-force scan of a few variables, for simpler cases. It is however not easy to derive more than about ten relations for the three difficult n 0 f parts of P (2) gq . For these cases we have employed the program provided in Ref.
[39], see also Refs. [40,41] to solve the remaining system of linear Diophantine equations. Since this program looks for short vectors, it is best for our purposes to eliminate 4 to 6 'unpleasant' coefficients, in particular those of low-weight combinations such as D 2 0 , D 2 1 , D 2 0 S 1 , D 2 1 S 1 , using the moments to N = 9 or N = 13, and work with the remaining 6 to 8 equations. For example, using the moments (3.1) this procedure leads to the result in the standard (M) definition of the MS scheme [56], where we have again used the abbreviations (3.4) and (3.6) and suppressed the argument N of the harmonic sums. The corresponding expressions for the C A C 2 F and C 2 A C F parts are somewhat longer, see below. The n f -dependent terms are much shorter; their determination does not require the N = 23 and N = 25 moments.
Note the simplicity of the coefficients in Eq. (3.8), in particular those of the terms with overall weights of 5 and 4 and sums of weight 2 or higher, which strongly indicates that the result is correct even without further checks. In fact, if any erroneous information is entered for an externally fixed parameters, e.g., a wrong coefficient of D 5 1 , or if the set of functions is too small, e.g., by omitting the term with D −1 , then either no solution exists for the system of Diophantine equations, or only solutions with nonsensically large coefficients (also) for the high-weight terms.
Nevertheless it is, of course, necessary to validate the resulting all-N formulae. For this purpose their predictions at higher values of N have been compared to additional MINCER moments such as The diagram calculations for the corresponding result at N = 29 have been carried out only in the planar limit C A − 2C F → 0 at n f = 0. As this result combines the three difficult all-N expressions for the C 3 F , C A C 2 F and C 2 A C F colour factors, which have been obtained independently from each other, it provides another strong check of all these results including Eq. (3.8). Perfect agreement is found for the not entirely trivial fractions at both values of N .
The overall most difficult case was the n f -independent, i.e., C 3 A part of P (2) gg . Also here the harmonic sums beyond weight 3 can be determined from the unpolarized case; the same holds for all terms not involving any un-summed denominators: these contribute to either the 1/(1 − x) + of the δ(1 − x) terms the large-x limit which are the same for P gg and P gg . This reduces the problem to the same basis set as in the case of P (2) gq at n f = 0. The 1/N 2 suppression of δ gg with respect to P gg , however, only removes one instead two coefficients for each harmonic sum up to weight 3.
Taking into account our additional knowledge of the coefficients of D 5 0 from Ref. [64] (this coefficient is the same for K (2) 66 and P (2) gg , unlike for the off-diagonal cases), D 5 1 and D 4 1 , cf. Refs. [68,69], and of S 1,1,1 , cf. Ref. [63], this leaves 49 terms with D 0 and D 1 plus the functions with the 'extra' denominator D −1 D 2 corresponding to D −1 in the previous case of P (2) gq .
The non-C F parts of P gg are non-singlet-like quantities, e.g., they are not affected by scheme transformations with z gg = 0, see Eqs. (2.18) and (2.19). Hence we could use some non-singlet heuristics, see Ref. [38], to reduce the overall basis to 52 functions, which we were able to determine using our own programs and, in the final step, Ref.
Quite a few of the resulting coefficients are far less simple than those in Eq. (3.8), see Eq. (4.12) below; on the other hand seven coefficients put in are zero, and there are some ex-pected relations. The result has been checked against the MINCER calculations at N = 27 and N = 29 which were finished only after we had obtained P (2) gg (N ). Another important check is the first moment which is not accessible directly [28], but can be obtained by Mellin-inverting to x-space expressions in terms of harmonic polylogarithms [71] from which arbitrary moments can be calculated. The results is see Eq. (2.16), as expected from the two previous orders. This result is the same in all factorization schemes considered here also for the C F terms due to P (n)

The NNLO splitting functions in Mellin space
The analytical odd-N expressions of the splitting functions to NNLO can be written in terms of harmonic sums [61] as recalled in Eqs. (3.2) and (3.3) above. Our notation is different from Section 3 of Refs. [12,13]: here all sums are taken at argument N (which we usually suppress), for the additional un-summed denominators we employ the abbreviations (3.4), (3.6) and In this notation the leading-order (LO) contributions [1,6,7] to Eq. (2.7), see also Eq. (2.9), read and their next-to-leading order (NLO) counterparts of Refs. [23][24][25] are given by For completeness also including the non-singlet contribution, which is identical to the function P −(2) ns (N ) given (in a different notation) already in Eq. (3.8) of Ref. [12], the polarized next-to-next-to-leading (NNLO) quark-quark splitting function P (2) qq (N ) is the sum of and P (2) ps In N -space the off-diagonal NNLO entries of the matrix (2.7) are given by Finally the polarized third-order gluon-gluon splitting function reads All these results refer to the standard transformation to the MS scheme of Ref. [56], see Eq. (2.17). With the exception of the C A n 2 f part of Eq. (4.12), which was derived in Ref. [73] (see also Ref. [74]), Eqs. (4.9)-(4.12) are new results of the present article.
The last two equations include the denominator ν defined in Eq. (4.1), and are therefore only valid at N ≥ 3. The first moment of the NNLO quark-gluon splitting function is 14) The corresponding results for the gluon-gluon splitting function are identical to the coefficients of the beta function recalled in Eq. (2.16). The NLO and NNLO pure-singlet results are related to Eqs. (4.13) and (4.14) by In the OPE, this relation for the anomalous dimension of the pure-singlet axial current together with Eq. (3.10) for the first moment of P (2) gg has been shown in Ref. [49] to be a direct consequence of the requirement that the axial anomaly [82,83] should preserve the one-loop character of the operator relation [84] ∂ μ j 5 μ = −2n f a s G μν a G a,μν (4.17) in dimensional regularization, where j 5 μ = ψγ 5 γ μ ψ and G μν a ( G μν a = 1/2ε μναβ G a,αβ ) denote the renormalized axial current and the (dual) gluon field-strength tensor. In this context Eqs. (3.10) and (4.16) are thus consistency requirements ensuring the correct renormalization of the pure-singlet axial current with the chosen finite renormalization constants Z ik , see Eq. (2.17). Consequently Eq. (4.16) for n = 3, together with Eq. (4.15) and P +(n) ns (N = 1) = P (n) qg (N = 1) = 0, fixes the first moments of the upper-row splitting functions at order α 4 s . The quantities given above do not provide the complete set of third-order helicity-difference splitting functions. Additional even-N functions P −,v ns exist for the quark-antiquark differences that occur in the (so far practically irrelevant) structure functions g 3 and g 4 in polarized chargedcurrent DIS which has been analyzed at NLO in Ref. [75]. The corresponding NNLO corrections may be addressed in a future publication together with the generalization of Refs. [76,77] to all N . It appears safe to assume P −(2) ns = P +(2) ns as given in Eq. (3.7) of Ref. [12], P v(2) ns is unknown though at this point.
The small-x behaviour in the right parts of Figs. 5-8 is due to the above contributions, which exhibit the usual pattern of alternating LL, NLL etc terms with coefficients strongly increasing towards lower logarithms. Consequently the leading logarithms alone do not provide a good approximation for any practically relevant values of x as illustrated in the figures. Yet it is also clear, from the scale of the ordinates in those right panels and Eq. (5.20)-(5.23), that these logarithms lead to a huge small-x enhancement that can potentially spoil the stability of the expansion in α s at x-values that would be accessible to an electron-proton collider with polarized beams.
Given the length and complexity of the exact expressions (5.12)-(5.16), it may be useful to also have at one's disposal compact and accurate approximate expressions for the case of QCD, i.e., C A = 3 and C F = 4/3. Such approximations can be build up, besides powers of x, from the non-logarithmic plus distribution and end-point logarithms ns , the result (4.23) of Ref. [12] can be used also here; it is given by  The other terms at x < 1 have been fitted to the exact results, evaluated by the FORTRAN code of Ref. [79], at 10 −6 ≤ x ≤ 1 − 10 −6 using the MINUIT package [80,81]. Except for x-values very close to zeros of the splitting functions, the above parametrizations deviate from the exact results by less than one part in thousand, which should be sufficient for any foreseeable phenomenological application. As in the unpolarized case [12,13], the coefficients of δ(1 − x) have been adjusted in Eq. (5.29) using low integer moments in order to achieve a maximal accuracy of the parametrization and its convolutions with the polarized gluon distribution. For a brief discussion of this slightly subtle point the reader is referred to Ref. [13] (penultimate paragraph of Section 4).
The effect of the new results (5.13)-(5.16) on the evolution of polarized parton densities is briefly illustrated in Figs. 9 and 10, where the respective first and second lines of Eq. (2.7) have been evaluated for the schematic, but sufficiently realistic low-scale distributions used for the evolution benchmarks in Refs. [85,86], for α s (μ 2 0 ) = 0.3 and n f = 3. After the convolution with the distributions (5.30), the NNLO corrections are fairly small down to small x.

Summary
We have extended the determination of the helicity-difference (polarized) splitting functions P ik , which were only known at the first [1,6,7] and second [23][24][25] order in the strong coupling constant α s so far, to the third order (next-to-next-to-leading order, NNLO) in massless perturbative QCD. These corrections are relevant to the structure function g 1 in polarized deepinelastic scattering (DIS), for which we also confirm the results of Ref. [44] for the NNLO coefficient functions, and all other observables that are sensitive to the polarized quark and gluon distributions f q i + fq i and f g . The so far practically irrelevant polarized quark-antiquark differences have not been addressed here; the corresponding splitting functions can be calculated, e.g., by extending the analysis of weak-interaction structure functions in Ref. [75] to NNLO accuracy.
The calculation of the upper row of the matrix of NNLO flavour-singlet splitting functions, i.e., of P (2) qq (x) and P (2) qg (x), was carried out via the structure function g 1 as a direct extension of our previous calculations of the helicity-averaged (unpolarized) case [12,13], for an earlier brief account see Ref. [26]. The corresponding lower-row quantities P (2) gq (x) and P (2) gg (x) have been determined in a different manner from graviton-exchange DIS, see Ref. [28], which includes structure functions sensitive to the polarized gluon distribution at the Born level.
We have first calculated the relevant structure function at fixed odd moments to N = 25, using a large-N optimized version [35] of the MINCER program [33,34] in (T)FORM [31,32]. Exploiting in particular the close relation between the polarized and unpolarized splitting functions for the highest-weight harmonic sums [61] and for the threshold limit, cf. Ref. [58] -which includes the so-called supersymmetric relation, see Refs. [3,87], as far as it can be addressed in MS -we have then been able to determine the all-N expressions of P (2) gq and P (2) gg . It was crucial for this step that the coefficients of the harmonic sums are integer, up to low powers of 2 and 3 that can be removed by a suitable normalization, which allows the use of advanced tools [39-41] for systems of Diophantine equations; this was observed and exploited before in a comparable but somewhat simpler situation in Ref. [38]. Finally the results have been validated by comparing the next two moments of all-N expressions to additional results calculated using MINCER up to N = 29.
Our results have been presented above in N -space and x-space, using the transformation of Ref. [56] from the so-called Larin scheme for γ 5 [48,49] in dimensional regularization to MS. This scheme shows an unphysical feature in the threshold limit of the quark-gluon splitting function P gq already at NLO, which can be removed to NNLO by simple additional terms in the scheme transformation. Yet this situation does not appear to necessitate a change of the factorization scheme in practical calculations after almost two decades of NLO analyses in QCD spin physics.
The new functions P (2) ik (x) are consistent with all known limits and partial results, e.g., for the leading large-n f terms [73], and expectations. In particular, the first moment of P (2) gg (x), which is not directly accessible in graviton-exchange DIS [28] but can be determined from the x-space results in terms of harmonic polylogarithms [71], is identical to the NNLO coefficient of the beta function of QCD [54,55] as theoretically required. We have checked our calculations of graviton-exchange DIS also by re-calculating, and obtaining full agreement for, P (2) qq and P (2) qg to fairly high values odd of N and all unpolarized flavour-singlet NNLO splitting functions at even N ≤ 10. As those results, the present polarized splitting functions lead to fairly small NNLO corrections, down to low values of x, after the convolution with realistic polarized quark and gluon distributions, despite a double-logarithmic small-x enhancement that dwarfs that of the non-singlet cases.
Our results allow NNLO analyses of spin-dependent hard-scattering observables, provided that the corresponding coefficient functions are known to this accuracy as for the structure function g 1 in DIS [44], for a fixed number of effectively massless flavours n f . The extension to analyses in the so-called variable flavour-number scheme, where effective theories for different values of n f are used together, requires non-trivial matching coefficients for the strong coupling [88] and the parton densities at this order. The latter coefficients have been calculated in Ref. [89] for the unpolarized case. As far as we know, the corresponding results for the helicity-difference parton distributions are not yet available in the literature though.
FORM and FORTRAN files of our main analytical results in N -space and x-space, and compact high-accuracy parametrizations of the functions P (2) ik (x), can be obtained by downloading the source of this article from http://arxiv.org/ or from the authors upon request.

Appendix C. Calculation of graviton-exchange DIS
Here we present some core ingredients of our diagram calculations, starting with the Feynman rules as used for graviton-exchange DIS. They have been taken from various sources [90,91]. We assume all momenta of the gluons and the graviton to be outgoing, while the momenta of the quarks and ghosts follow the arrows on the lines. The color indices in the fundamental representation are i and j ; color indices in the adjoint representation are represented by the letters a, b, c, d, e; the Lorentz indices of the graviton are α and β and those of the gluons are μ, ν, ρ, σ . We also use a gauge parameter which is indicated by ξ .
In the polarized case we do not need a ghost contribution, neither for the graviton nor for the gluon. Propagators for the graviton and the corresponding ghost are not required since we do not consider internal gravitons.