Four-Loop Non-Singlet Splitting Functions in the Planar Limit and Beyond

We present the next-to-next-to-next-to-leading order (N^3LO) contributions to the non-singlet splitting functions for both parton distribution and fragmentation functions in perturbative QCD. The exact expressions are derived for the terms contributing in the limit of a large number of colours. For the remaining contributions, approximations are provided that are sufficient for all collider-physics applications. From their threshold limits we derive analytical and high-accuracy numerical results, respectively, for all contributions to the four-loop cusp anomalous dimension for quarks, including the terms proportional to quartic Casimir operators. We briefly illustrate the numerical size of the four-loop corrections, and the remarkable renormalization-scale stability of the N^3LO results, for the evolution of the non-singlet parton distribution and the fragmentation functions. Our results appear to provide a first point of contact of four-loop QCD calculations and the so-called wrapping corrections to anomalous dimensions in N=4 super Yang-Mills theory.


Introduction
Within the gauge theory of the strong interaction, Quantum Chromodynamics (QCD), the precision of theory predictions for hard reactions at colliders crucially depends on our knowledge of hadronic matrix elements for the description of the long-distance hadronic degrees of freedom, once the hard-interaction part due to short-distance physics has been separated by means of QCD factorization. For scattering reactions with initial-state protons the relevant matrix elements are given by the well-known parton distribution functions (PDFs) of the proton, which provide information about the fractions of the proton's longitudinal momentum carried by the partons.
The dependence of these PDFs on the scale Q 2 is generated by evolution equations for the corresponding local operator matrix elements (OMEs). The relevant anomalous dimensions as functions of the Mellin moment N , or splitting functions as functions of the momentum fraction x, can be computed order by order in perturbative QCD. The corresponding one-and two-loop results have been known since long [1][2][3][4][5][6][7][8][9][10][11][12][13]. The current precision is at the three-loop level [14,15] -see refs. [16][17][18][19] for partial recalculations of these results -i.e., at the next-to-next-to-leading order (NNLO), which is nowadays the accepted standard for analyses of PDFs [20] and forms the backbone of precision predictions at the Large Hadron Collider (LHC).
However, computations for a number key observables at hadron colliders have been performed even at next-to-next-to-next-to-leading order (N 3 LO), including the cross section for Higgs-boson production in gluon-gluon fusion [21] and structure functions in deepinelastic scattering (DIS) [22][23][24][25]. The latter results have also found an application in predicting Higgs-boson production in vector-boson fusion at the LHC [26]. Due to QCD factorization, the resulting predictions carry a residual uncertainty and dependence on the factorization scheme due to the missing N 3 LO (i.e., four-loop) splitting functions. This situation motivates the computation of the QCD splitting functions at four loops. First steps in this direction have already been taken in refs. [27][28][29][30][31] at low N , and in ref. [32] where large-n f contributions have been derived at all N .
In the present article, we address the splitting functions for the non-singlet quark evolution equations at four loops in QCD. We use Forcer [33], a Form [34][35][36] program for four-loop massless propagators, to compute the anomalous dimensions at fixed integer values of the Mellin variable N . In the planar limit, i.e., for large n c for a general colour SU(n c ) gauge group, the exact four-loop results for moments up to N = 20 turn out to be sufficient to find and validate the analytic expressions as functions of N in terms of harmonic sums [37,38] by LLL-based techniques [39,42] 1 for solving systems of Diophantine equations. Such an approach has been used for anomalous dimensions at the three-and four-loop level before, cf. refs. [32,43,44]. Our analytic results in the threshold limit x → 1 ( N → ∞) include the (light-like) four-loop cusp anomalous dimension, see ref. [45], which has also been obtained in refs. [46,47] by different means.
Beyond the large-n c limit, we have computed the moments up to N = 16 for a general gauge group. These results are insufficient for a reconstruction of the analytic all-N results. They can be used, though, to obtain approximations for the four-loop splitting functions including x-dependent estimates of their residual uncertainties, see, e.g., earlier work at the three-loop level [48][49][50]. The approximations presented below are sufficiently accurate for the evolution of non-singlet PDFs down to small x, and include numerical results for the non-planar contributions to the four-loop cusp anomalous dimension that are sufficiently precise for phenomenological applications.
For processes with identified hadrons in the final state, QCD factorization requires fragmentation functions (FFs) that account for the physics of hadronization at long distances. Completely analogous to PDFs, the scale dependence of FFs can be computed within perturbative QCD. However, in contrast to the case of initial state hadrons, where the evolution equations for the scale-dependence of the PDFs are controlled by space-like kinematics, Q 2 ≤ 0, the scale evolution of the FFs with Q 2 ≥ 0 requires the so-called timelike splitting functions. These functions are known completely at two loops [9][10][11][51][52][53], see also refs. [54,55]. The three-loop corrections have been obtained in refs. [56][57][58] up to a phenomenologically irrelevant small uncertainty in the result for the time-like NNLO quarkgluon splitting function. First NNLO analyses of FFs have been performed recently [59,60].
The three-loop results in refs. [56][57][58] have been derived using well-known relations between space-and time-like kinematics, i.e., the Drell-Yan-Levy relation for the analytic continuation in energy q 2 → −q 2 and the Gribov-Lipatov relation in x-space [61,62], see also refs. [63,64], and generalizations based on conformal symmetry yielding a universal reciprocity-respecting evolution kernel [65][66][67]. Exploiting these relations, it is possible to use (space-like) DIS results to predict (time-like) cross sections for single-particle inclusive electron-positron annihilation. Thus, we are able to present here also the flavour non-singlet evolution equations for FFs at four loops in QCD.
This article is organized as follows. In section 2 we specify our notations and present the theoretical framework for obtaining our results. In particular we address the basis of non-singlet operators, their renormalization and the respective anomalous dimensions. We sketch the work-flow of the perturbative computation up to four loops, list all colour factors to this order and discuss general and end-point properties of the anomalous dimensions and splitting functions.
In section 3 we present the results of our fixed-N diagram calculations of the fourloop non-singlet anomalous dimensions and their all-N generalization in the large-n c limit. We discuss the large-N behaviour of the latter which includes the four-loop cusp anomalous dimension. The x-space counterparts of these anomalous dimensions, i.e., the splitting functions, are addressed in section 4. We present the exact formulae and compact parametrizations for the large-n c splitting functions, and approximate expressions for all cases that cannot be obtained exactly for now.
Two important applications of these results are presented in section 5: we present high-accuracy numerical results for large-x coefficients, in particular the four-loop cusp anomalous dimension in QCD, and illustrate the N 3 LO evolution of all three types of nonsinglet quark distributions. The N 3 LO non-singlet evolution is extended to the 'time-like' case of final-state fragmentation functions in section 6. We summarize our main results and provide a brief outlook in section 7.
The appendices contain the Feynman rules in appendix A, the exact results for the anomalous dimensions at 1 ≤ N ≤ 16 at four loops in appendix B, and the analytic expression for the difference of the time-like and space-like four-loop splitting functions in appendix C. Finally appendix D provides the complete all-N result for the terms with ζ 5 , which may be of theoretical interest.

Theoretical framework and calculations
The standard set of spin-N twist-two irreducible flavour non-singlet quark operators is given by where ψ represents the quark field, D µ = ∂ µ −igA µ the covariant derivative, and λ α the diagonal generators of the flavour group SU(n f ). It is understood in eq. (2.1) that the symmetric and traceless part is taken with respect to the Lorentz indices µ i in the curly brackets. We consider (spin-averaged) matrix elements of these operators (OMEs), specifically for external quark (or anti-quark) fields with momenta p 1 and p 2 . The operators O ns in eq. (2.2) are contracted with tensors of rank N , where ∆ is a light-like vector, ∆ 2 = 0. In the present context we need to compute OMEs of renormalized operators with zero momentum flow through the operator vertex, thus p 1 = p 2 = p in eq. (2.2) for the (off-shell, p 2 = 0) momenta of the external (anti-) quarks, Here and below we use square brackets [. . . ] to denote renormalized operators (in a minimal subtraction scheme [68,69] of dimensional regularization [70,71]). This reduces the vertex diagrams for the OMEs to quark two-point functions and, therefore, the computational complexity to propagator-type diagrams. The perturbative expansion of the operator in eq. (2.1) contracted with eq. (2.3) generates vertices with additional gluons as depicted in figure 1. The current four-loop calculation requires up to four additional gluons. The corresponding Feynman rules are presented in appendix A, see refs. [5,72] for earlier calculations at two-and three-loop accuracy. In order to derive the anomalous dimensions for the scale dependence of the non-singlet PDFs we need to perform the renormalization of those operators O ns , which proceeds multiplicatively as The anomalous dimensions γ ns governing the scale dependence of these operators, are connected to the factors Z ns in eq. (2.5) by All flavour differences of quark-anti-quark sums (+) and differences (−) evolve in with the same anomalous dimensions γ + ns (N ) and γ − ns (N ), and the total valence distribution with γ v ns (N ) = γ − ns (N ) + γ s ns (N ), see, e.g., ref. [14]. These quantities are related to the corresponding splitting functions P ± ns (x) and P s ns (x) by a Mellin transform, where the relative sign is a standard convention. In perturbation theory these quantities can be expanded in powers of the strong coupling constant α s . Here and below we normalize a s ≡ α s /(4π), so that up to four loops 9) and similarly for the splitting functions P ns (x) and other quantities. The first-order quantity γ (0) ns is the same for all three cases given above. γ + ns and γ − ns differ at order a 2 s , and a non-vanishing flavour-independent ('sea') contribution γ s ns occurs at order a 3 s for the first time [14]. The fourth-order contributions γ (3) ns (N ) to all three quantities are addressed in the present article.
The actual computation follows a well-established production chain. The Feynman diagrams for the OMEs in eq. (2.2) are generated up to four loops using Qgraf [73]. The latest version [74] of the symbolic manipulation program Form [34,35] and its multithreaded version TForm [36] are used for all further steps. The Qgraf output is processed by a program that assigns the topology and computes the colour factor using the code of ref. [75]; the group invariants occurring in the present case are listed in table 1. Diagrams of the same topology and colour factor are combined to meta diagrams for computational efficiency, where lower-order self-energy insertions are treated as described in ref. [76]. Considering all color factors, this procedure leads to 1 one-loop, 7 two-loop, 53 three-loop and 650 four-loop meta diagrams for γ ± ns ; and 1 three-loop and 29 four-loop meta diagrams for γ s ns . For comparison: the output of Qgraf consists of 15901 four-loop diagrams. The running of the meta diagrams is managed using the database program Minos [77].
The diagram calculations are done in dimensional regularization [70,71] with the Forcer program [33] which was already used for the N ≤ 6 and high-n f computations in refs. [31,32]. Our agreement (after renormalization, see below) with those results, which were obtained in a different theoretical framework, provides a strong check of our present setup. The Forcer program itself has been validated in calculations of the four-loop renormalization of Yang-Mills theories to all powers of the gauge parameter, see ref. [79], and has recently been applied -together with the algorithms for the R * operation [80,81] developed in ref. [82] -in five-loop computations of the beta function, Higgs-boson decays to hadrons and the R-ratio in e + e − -annihilation in refs. [83,84].
The bare results A ns for the OMEs in eq. (2.4) obtained in this way are then subject to renormalization which we perform in the standard modified minimal subtraction scheme MS [68,69]. In this scheme and in D = 4 − 2ǫ dimensions the strong coupling α s evolves according to  Table 1. The colour factors for non-singlet OMEs up to four loops with their numerical values in SU(n c ) with N R = n c and QCD, see also ref. [78] for a discussion on the normalization of d abc d abc . The second column gives the notations of the result files, which are distributed with this article on https://arxiv.org. As in many other articles, we suppress the colour factor T F ( = 1/2 in SU(n c )) which can be readily re-instated.
where β(a s ) denotes the usual four-dimensional beta function in QCD, with coefficients β 0 = 11/3 C A − 2/3 n f etc, and n f represents the number of active quark flavours.
Using eq. (2.5) the renormalized OMEs [A ns ] are obtained by where we have made all dependences on N explicit. The factor Z ψ denotes the quark wave function renormalization constant accounting for the external quarks field with off-shell momenta in eq. (2.4), see, for instance ref. [79]. Unlike Z ns , the quantities Z ψ and A ns are gauge-dependent, hence also the renormalization constant Z ξ of the gauge parameter is required. The resulting operator renormalization factors Z ns in eq. (2.5) can be expressed as a Laurent series in ε as ns . (2.12) In this manner, the anomalous dimensions γ ns have been computed for a general gauge group at 1 ≤ N ≤ 16, i.e, γ + ns at even N and γ −,s ns at odd N . The exact results are listed in appendix B; numerical values for QCD can be found in section 3. The hardest (nonplanar) diagrams do not contribute in the limit of a large number of colours n c , where the functions γ + ns (N ) and γ − ns (N ) are identical, as it is evident from diagrammatical analyses and the known x-space expressions for P ± ns (x), see refs. [14,32,66,85]. Consequently we were able to obtain the even-N and odd-N values of the large-n c anomalous dimension, which is structurally simpler than full QCD results, even up to N = 20.
So far, fixed-N values of anomalous dimensions have been found to be fractions of (large) integer numbers, multiplied at most by values ζ 3 . . . ζ 2L−3 of the Riemann zetafunction at L loops. The denominator structure of the fractions suggests analytic all-N expressions in terms of harmonic sums [37,38] up to weight 2L − 1. Assuming no numerator-N terms, cf. refs. [22][23][24], the most complicated parts (without a factor ζ n ) of the non-singlet anomalous dimensions at n loops read where D k a are simple denominators, and S w (N ) is a shorthand for all harmonic sums of a given weight w with S 0 (N ) ≡ 1.
The calculated moments suggest a = 0, 1 for γ (3)± ns (N ), as at three loops [14] and for the n 2 f and n 3 f four-loop contributions [32]. The function γ s ns (N ), on the other hand, includes terms with a = −1 and a = 2.
The functions γ (3) ns (N ) contain harmonic sums up to weight w = 7, hence the ansatz (2.13) includes far too many unknown coefficients c akw for a direct determination from the (small) number of calculated moments. However, these coefficients are integer modulo some predictable powers of 2 and 3. Therefore the systems of equations derived from eq. (2.13) can be turned into Diophantine systems which require far fewer equations than unknowns and which can be solved by LLL-based techniques [39][40][41][42]. This approach has been successfully applied before in refs. [32,43,44].
In this context it is crucial to constrain eq. (2.13) as far as possible based on general properties of the anomalous dimensions γ ns (N ). Here three issues are worth pointing out. First, the functional forms of the γ ns (N ) are (conjectured to be) constrained by 'self-tuning' [66,67], where σ = −1 (+1) for the space-like (time-like) anomalous dimensions, and the nonsinglet universal evolution kernel γ u is reciprocity-respecting (RR), i.e., invariant under the replacement N → (1 − N ). By expanding the r.h.s. of eq. (2.15) about N and inserting the perturbation series of all quantities involved, γ u can be expressed in terms of the MS anomalous dimensions, see also ref. [86]. Expressing the latter in terms of γ 0 = γ S(n) ns = γ T (n) ns and the average of the space-like and time-like expansion coefficients γ n = , one arrives at where we have used the abbreviation d N = d/dN and suppressed the N -dependences for brevity. A convenient way to take these derivatives is via inverse Mellin transforms to x-space, where the multiplication with ln n x corresponds to the N -space operator d n /dN n , and Mellin transforms of the result. The required manipulations can be readily performed using algorithms for harmonic sums, harmonic polylogarithms and their (inverse) Mellin transformations [37,87,88] which have been implemented in publicly available Form packages described in ref. [34].
Since the difference between the time-like and space-like anomalous dimensions is known to four loops, eq. (2.13) can be applied to the RR quantity γ ns . This implies that the denominators 1/N and 1/(N + 1) can only enter in the combination 1/(N (N + 1)), and that only RR (combinations of) harmonic sums occur, see refs. [89,90], which reduces the number of sums at weight w from 2 · 3 w−1 to 2 w−1 . Assuming that only powers of 1/(N + 1) enter in addition, the total number of basis functions in eq. (2.13) up to weight w is 2 w+1 − 1, e.g., 255 for w = 7. Even taking account end-point constraints, see below, this is a prohibitively large number for now.
Second, the identical leading-n c terms of γ ± ns (N ) contain only non-alternating harmonic sums, i.e., only positive indices in eq. (3.4). This reduces the number of RR sums of weight w to the Fibonacci number F (w), i.e., 1, 1, 2, 3, 5, 8, 13 for w = 1 to w = 7, as can be seen by counting the number of binomial harmonic sums at weight w [89]. Considering all combinations with additional powers of the weight-1 object 1/(N (N + 1), the total number of functions up to weight w in eq. (2.13) amounts to F (w + 4) − 2, e.g., 87 for w = 7.
The third and final point is that the N → ∞ (large-x) and N → 0 (small-x) limits of the anomalous dimensions (splitting functions) provide a substantial number of constraints. If one disregards terms of order O(1/N 2 ) for N → ∞ , then all three non-singlet anomalous dimensions γ a ns (N ), a = +, −, v, are identical and given by [65] (γ e is the Euler-Mascheroni constant) Here the coefficients A n -the n-loop (light-like) cusp anomalous dimension -and B n provide genuine n-loop information. The coefficients C n and D n , on the other hand, can be expressed in terms of lower-order information (see eqs. (3.10) and (3.11) below). This and the absence of second and higher powers of ln N in eq. (2.17), and similar if less stringent constraints on N −k ln N terms with k > 1, provide a substantial number of constraints on the coefficients in eq. (2.13).

JHEP10(2017)041
The small-x expansion of the splitting function P (n) ns (x) shows a double-logarithmic enhancement, i.e., there are contributions of the form x a ln b x with a > 0 and b ≤ 2n. The leading-logarithmic (LL) contributions to P ± ns have been known to all orders for a long time [91,92]. This resummation has been extended to next-to-next-to-logarithmic accuracy for the x 2k ln b x contributions to P + ns and the x 2k+1 ln b x contributions to P − ns at all k ≥ 0 [93,94]. The formal structure of these results is analogous to their time-like counterparts [95,96], but the numerical pattern is completely different such that the spacelike resummation is of no direct phenomenological use. The functions P + ns (x) and P − ns (x) are the same in the large-n c limit, hence in this case the small-x resummation constrains the coefficients contributing to for a ≥ 0 and 4 ≤ b ≤ 6 .
(2.18) An alternative approach to the limit N → 0, i.e., the small-x logarithms for a = 0, has been pursued in ref. [97]. In the large-n c limit, the generalization of the LL relation in refs. [91,92] correctly (re-) produces all a = 0 small-x logarithms obtained in refs. [14,32,93], after correcting typos in eqs. (25) and (26) of ref. [97]. Hence we can assume that eq. (2.19) is also correct for the n 0 f and n 1 f four-loop contributions at large n c .
Together these relations comprise 18 large-N and 28 small-x constraints for the n 0 fterms at four loops eliminating more than half of the 87 free parameters of the w = 7 large-n c ansatz, after which it is possible to solve the remaining system of Diophantine equation using the moments N = 1, . . . , 18 with the program axb() of the Calc package [42]. The resulting analytic expressions for γ (3) ns (N ) agree with the result of the diagram calculations at N = 19 and N = 20. This agreement renders it extremely likely -although, of course, not mathematically certain -that these results (and, therefore, the above structural conjectures and features used in their derivation) are correct.
As mentioned above, present information and understanding appears not to be sufficient for extending these analytic results beyond the large-n c contributions for the n 0 For the remaining functions we resort to x-space approximations based on the first eight even-N or odd-N moments supplemented by the large-x and small-x constraints discussed above. These approximations and their error estimates can be constructed in the same manner as those for the three-loop splitting functions in refs. [48][49][50]. The present results are more accurate, though, due to the higher number of available moments and the improved understanding of the end-point limits. The fact that the large-N limit (2.17) includes only the two free parameters A 4 and B 4 , in particular, results in a high accuracy of these coefficients which are relevant also in the context of the soft-gluon exponentiation, see refs. [98][99][100][101] and references therein, and beyond.

Results in N -space
We start presenting our results by writing down the moments to N = 16 of the non-singlet four-loop anomalous dimensions for QCD in a numerical form. The exact results for a general gauge group with one set of fermions can be found in appendix B. For γ (3)± ns (N ) we separately display the leading (subscript L) and non-leading (subscript N ) contributions in the large-n c limit of SU(n c ) at n c = 3. The former correspond to the colour factors C F n 3−k c n k f . The latter collect all other terms, which are suppressed by two or more powers of n c , cf. table 1 above.
The first eight even-N values of γ It is clear from these results, that the large-n c limit alone provides an excellent approximation to the individual n a f coefficients except for the lowest values of N . The non-large-n c 'correction' amounts to 10% and 4% for the n 0 f and n 1 f terms, respectively, but 2% or less at N ≥ 7 in both cases.
We have computed the first nine odd-N values of the 'sea' contribution γ in the large-n c limit. The complete lowerorder contributions can be found, in a different notation but the same normalization, in eqs. (3.4) -(3.8) of ref. [14]. The anomalous dimensions can be expressed in terms of the denominators D k a in eq. (2.14) and harmonic sums [37,38] at argument N , which are recursively defined by The weight w of the harmonic sums is defined by the sum of the absolute values of the indices m d . Sums up to w = 2 n − 1 occur in the n-loop anomalous dimensions. The argument N of the sums is suppressed for brevity below, and we use the shorthand η = 1/(N (N + 1)) = D 0 D 1 . (3.7) The large-N limit of eq. (3.5) is of the form (2.17) with the large-n c cusp anomalous dimension Our result for the (complete) n 2 f part was first presented at Loops & Legs 2016, see ref. [31], the rest in a Zurich seminar by one of us [103]. Eq. (3.8) agrees with results of refs. [46,47], where this quantity was obtained by computing the photon-quark form factor in the largen c limit. The lower-order coefficients can be found in eq. (3.11) of ref. [14].
The one-to three-loop coefficients B 1,2,3 in eq. (4.9) can be found, as coefficients of δ(1−x), in eqs. (4.5), (4.6) and (4.9) of ref. [14]. The four-loop coefficient in the large-n c limit reads The coefficients B n contain collinear contributions to the evolution kernels. With the help of the QCD corrections to the quark form factor in dimensional regularization, one can extract from them the universal eikonal anomalous dimension. The latter governs the subleading infrared poles in gauge-theory amplitudes and captures contributions from large-angle soft gluons [100,104,105]. As mentioned above, the coefficients C n and D n in eq. (4.9) do not provide new information, but are functions of lower-order quantities. They are given by cf. ref. [65], which leads to the four-loop relations Using the results (3.8) and (3.9), it is actually now possible to predict C 5 and D 5 for large n c . The new functions (3.6) and (3.7) are shown in figures 2 and 3, respectively, together with their large-N approximation (2.17) with the coefficients given above. In the right panels, the results are divided by ln N , so for N → ∞ the curves tend to constants given by the respective terms in the four-loop cusp anomalous dimension (3.8).
The approach to this asymptotic behaviour is very slow: the n 0 f contribution in figure 2 is 0.856 of its asymptotic result at N = 30, yet it deviates by less than 10% only from an N -value above N = 500. The corresponding numbers for the n 1 f part in figure 3 are (N ) in the large-n c limit, compared with its large-N expansion with all terms included in eq. (2.17) and, in the right panel, its asymptotic behaviour for N → ∞. The exact curve has been computed by via the x-space counterpart of eq. (3.6), see section 4. 0.873 at N = 30 and N ≃ 185 for a deviation by less than 10%. It might be interesting to note, on the other hand, that the corresponding n a f coefficient of A n , here and in all lower-order cases (in full QCD), falls in the interval spanned by the corresponding results for γ  figure 4, together with the all-N results in the large-n c limit. As at the previous orders in α s , there are cancellations between the n findependent and the n f -dependent contributions, which are particularly pronounced here at n f = 5. For this number of light flavours, which is relevant for high-energy processes at the LHC, the large-n c result do not describe the (small) fourth-order QCD contributions to the non-singlet evolution equations at the phenomenologically most relevant moments N and momentum fractions x. We therefore need to convert the calculated moments to practically usable constraints on the four-loop splitting functions P

Results in x-space
The four-loop non-singlet splitting functions P    performed by an algebraic procedure [87,88] based on the fact that harmonic sums occur as coefficients of the Taylor expansion of HPLs.
For the convenience of the reader, we recall their basic definitions [87]. The lowestweight (w = 1) functions H m (x) are given by The higher-weight (w ≥ 2) functions are recursively defined as For chains of indices 'zero' we employ the abbreviated notation The argument x will be suppressed in all results below, and we express the terms with (1 ± x) −1 in terms of the x-dependence of the leading-order splitting function P In this notation, the common large-n c limit of the functions P with the new results and The n 2 f and n 3 f contributions to eq. (4.6) are given by eqs. (4.6) and (4.12) of ref. [32]. Disregarding terms that vanish for x → 1, the large-x behaviour of P (3) ns, L (x) can be written as in terms of the coefficients specified in eqs. (3.8), (3.9) and (3.11) above. The numerical values of the coefficients of the small-x logarithms, ln k x with k = 1, . . . , 6, can be read off from eq. (4.11) below. All six logarithms and the constant contribution for x → 0 are required for a good approximation to the splitting functions at small x-values relevant for collider physics.
In view of the length and complexity of the exact expressions (4.7) and (4.8) it is useful to have at one's disposal also compact approximate representations involving, besides powers of x, only simple functions like the plus-distribution and the end-point logarithms (4.10) Such approximations can be readily used in N -space evolution programmes, see, e.g., ref. [106]. The results (4.7) and (4.8) can be parametrized with a high accuracy (of 0.1%

JHEP10(2017)041
or better) as Here the exact large-x and small-x coefficients have been rounded to seven significant figures. The brackets multiplied by 25000 and 2500 have been obtained by fits to the exact expressions at 10 −6 ≤ x ≤ 1 − 10 −6 . The small shifts of the δ(1 − x) fine-tune the accuracy of the resulting low moments and of the convolutions with the quark distributions. The required evaluation of the HPLs has been performed using a weight-6 extension of ref. [107] and the program of ref. [108], which return identical results at the accuracy considered here.
For the corresponding non-leading contributions in the large-n c limit, denoted by the subscript N in eqs. (3.1) and (3.2) above, we are for now limited to approximations analogous to (but more accurate than) those once constructed at three loops [48][49][50]. For the n 0 f and n 1 f parts of P • two of three suppressed large-x logarithms (1−x) ln k (1−x), k = 1, 2, 3, • one of ten two-parameter polynomials in x that vanish for x → 1, • two of the three unknown small-x logarithms ln k x, k = 1, 2, 3.
The parameters of the 90 resulting trial functions are determined from the eight available moments, and then two representatives as chosen that indicate the remaining uncertainty. The result of this process is illustrated in figures 5 and 6. Supplementing the approximations A and B in figures 5 and 6 by accurate parametrizations of the complete n 2 f results of ref. [32] and the exact (but numerically truncated) n 3 f expressions in a non-HPL notation, we obtain  and (4.14) The case of P (3)− ns (x) can be treated in the same manner, but taking into account that only its leading small-x logarithm is known up to now [92]. After careful consideration, the two approximations indicating the uncertainty band in this case are chosen as L,0 (x) + C F n 2 c n f P L,1 (x) + P L,1 (x) + P The n 3 f contribution to this last equation is the same as in eq. (4.14). Before we illustrate these results, it is useful to briefly recall the behaviour of the corresponding third-order splitting functions. This is done in figure 7 for n f = 4 flavours.
The corresponding size and uncertainty bands of P the corrections as coefficients of α 3 s and α 4 s , respectively, are comparable in the region of x for which definite conclusions can be drawn.
The four-loop 'sea' contribution P (3)s ns (x) to the evolution of the total valence distribution is suppressed by two powers of (1−x) for x → 1, but its n 1 f part is completely unknown in the small-x limit. In this case, we use the nine odd moments (3.3) with a suitably modified ansatz, in which the coefficient of ln 6 x is varied 'by hand' over a sufficiently wide range, and the coefficients of ln k x, k = 1, . . . , 5 , are all determined from the moments. In this manner we obtain The last equation is a high-accuracy parametrization, constructed in the same manner as eq. (4.11) above, of the exact result given in eq. (4.11) of ref. [32]. The trial functions considered for all three cases lead to very similar predictions for the respective next moments, i.e., N = 18 for P  ns (x) by deriving less accurate approximations using one fewer moment and comparing the results to the now unused highest calculated moments. As far as we can see from this and other checks, our approximation procedure, which is of course not mathematically rigorous, does not underestimate the remaining uncertainties.

Numerical implications
We are now ready to address two important applications of our new fourth-order results. First, as already mentioned above, the large-x / large-N limits of the splitting functions include coefficients that are relevant beyond the evolution of parton distributions: the (light-like) four-loop cusp anomalous dimension A 4 and the δ(1−x) coefficient B 4 for quark fields. We are now able to provide approximate if rather accurate numerical results for these coefficients. The obvious second application is a (further) improvement of the perturbative stability of the evolution of the non-singlet quark distributions over a wide range in x.
The analytic large-n c expression for A 4 has been presented in eq. (3.8) above. Together with the approximate results in eqs.  [98]. The agreement of the actual results with these approximants would be (much) better without the contributions of the quartic group invariant (see below). A similar situation has been observed for the four-loop beta function in ref. [109]. The expansion of the cusp anomalous dimension, now to the fourth order in α s , is given by the very benign series colour factor It may be interesting, for theoretical purposes, to consider the contributions of the individual colour factors to A 4 and B 4 . By repeating the approximation procedure of the previous sections separately for each colour factor, we arrive at the corresponding results collected in table 2. Our results show that both quartic group invariants definitely contribute to the four-loop cusp anomalous dimension -an issue that has attracted some interest, see, e.g., refs. [110][111][112][113][114] -which means that the so-called Casimir scaling between the quark and gluon cusp anomalous dimensions, A q = C F /C A A g , does not hold beyond three loops. A lower value, -113.66 after conversion to our notation, results from assumptions made in ref. [115] for the coefficient of n f d abcd We now turn to the effect of the four-loop splitting functions (4.6)-(4.21) on the evolution -specifically the logarithmic derivativesq i ns ≡ d ln q i ns /d ln µ 2 f where µ f is the factorization scale -of the non-singlet combinations q ±,v ns (x, µ 2 f ) of the quark and anti-quark distributions. In all three cases we employ the same schematic, but characteristic model distribution This facilitates a direct comparison of effects of the various contributions of the splitting functions. For the same reason the reference scale is specified by the order-independent value α s (µ 2 0 ) = 0.2 (5.7) for the strong coupling constant. This value corresponds to µ 2 0 ≃ 25 . . . 50 GeV 2 for α s (M 2 Z ) = 0.114 . . . 0.120 beyond the leading order. In this region of the physical scale Q 2 deep-inelastic scattering has been measured both at fixed-target experiments and, for much smaller x, at the ep collider HERA. Our default for the number of effectively massless flavours is n f = 4.
The reliability of perturbative calculations can be assessed by the relative size of the higher-order correction at a 'nominal' value of the renormalization scale µ r , here µ r = µ f , and by investigating the stability of the results under variations of µ r . For µ r = µ f the inverse Mellin transform, see eq. (2.8), of the perturbative expansion (2.9) in terms of a s = α s /(4π) has to be replaced by For the MS expansion coefficients β k of the beta function of QCD to N 3 LO see refs. [109,116] and references therein. In figure 9 the consequences of varying µ r over the range

JHEP10(2017)041
Outside the region around x = 0.07 where the µ f -derivatives are small, the remaining uncertainty ofq + ns is well below 1% down to x ≃ 10 −3 and possibly, below. The size and µ r -variation of the NLO and NNLO contributions are somewhat larger forq − ns at small x, yet neither the N 3 LO correction nor its scale variation exceeds 1% in the region shown in the plot.
The case of q v ns , shown in figure 12 is noticeably different beyond NLO due to the appearance of the d abc d abc n f contribution P s ns which is negligible and at large x but large at small x [14]: the difference of the NNLO curves of figure 11 and 12 -note the different scales for the ordinate -is caused completely by this contribution due to our choice (5.6) of the input quark distribution. Also in this case our new N 3 LO results leads to a considerable improvement and a remaining small-x uncertainty of no more than about 2% at x ≥ 10 −4 .

The time-like case
The differences between the initial-state ('space-like', σ = −1) and the final-state ('timelike', σ = 1) splitting functions respectively governing the evolution of the parton distributions and fragmentation functions can be expressed in terms of lower-order quantities. At N 3 LO they read where we have used the short-hand notations A ⊗2 ≡ A ⊗ A etc for the Mellin convolutions, and P (n) i stands for the average of the corresponding σ = 1 and σ = −1 expansion coefficients, normalized as in eq. (2.9). Eq. (6.1) has been derived in ref. [56] by generalizing results in ref. [65]; it is also a direct consequence of eq. (2.15) [67]. The resulting rather lengthy explicit expressions can be found in appendix C. Here we present parametrizations in terms of powers of x and the logarithms in eq. (4.10). As above, their small-x and large-x coefficients are exact up to their rounding to seven digits. The accuracy of the n k f coefficients is better than 0.1% except close to zeros. The three parametrizations are given by +xL 0 (13.301+2.1060L 0 +0.4375L 2 0 )+L 1 (13.060x 1 +11.023L 0 ))−9.550559L 0 The difference of the time-like and space-like splitting functions and the resulting scale derivative are illustrated in figures 13 and 14 for the most important case, NS + . The pattern is somewhat different in the time-like case, e.g., the N 2 LO contribution is negative at small x. Yet also here the perturbative expansion is 'perfectly' stable after including the N 3 LO corrections, with a residual uncertainty of 1% or less for x > 10 −4 at our (for the time-like case, low-scale) reference point.

Summary and outlook
We have presented the four-loop corrections P (3)±,v ns (x) to all three non-singlet splitting functions in perturbative QCD. Our results, which are partly approximate but sufficiently accurate for collider-physics applications, allow to set-up and solve the QCD evolution equations for flavour non-singlet (initial-state) parton distributions (PDFs) and (final-state) fragmentation functions (FFs) at N 3 LO. They thus provide a major step towards the consistent application of QCD factorization to theoretical predictions for N 3 LO cross sections with initial (final) state hadrons, as already obtained in refs. [21][22][23][24]26], which requires hard partonic cross section and PDFs (FFs) at the same order in renormalization-group improved perturbation theory.
The resulting logarithmic scale derivatives d ln q ±, v ns /d ln µ 2 f exhibit a very good convergence of the perturbative expansion. Both the four-loop corrections and the N 3 LO dependence on the renormalization scale µ r mostly amount to as little as 1% or less (and maximally 2%, for P v ns (x) at small x) at momentum fractions x > 10 −4 for α s (µ 2 f ) ≃ 0. 1 Figure 13. The perturbative expansion of the difference δ P + = P + ns, σ=1 − P + ns, σ=−1 of the timelike (σ = 1) and space-like (σ = −1) non-singlet + splitting functions. Left: results in Mellin-space, compared to beyond-LO contributions in the space-like case. Right: convolution with the schematic input shape (5.6).
For the large-n c contributions to P (3)± ns , these moments turn out to be sufficient for a reconstruction of the all-N results in terms of harmonic sums by solving systems of Diophantine equations. Additional knowledge about the limits for x → 0 and x → 1, as well as the rephrasing (based on conformal symmetry) of the evolution equations in terms of a universal 'reciprocity-respecting' evolution kernel [65][66][67], have been instrumental in this step. Beyond the large-n c contributions we have used the computed Mellin moments, again supplemented by endpoint constraints, to provide approximations for the four-loop splitting functions including x-dependent estimates of their residual uncertainties. The latter are very small, except in the region x < 10 −2 . These small-x uncertainties are subject to a further suppression in the actual evolution due to the convolution with the PDFs (or FFs). Due to this our results are found to be sufficiently precise for x 10 −4 .
From the threshold limit x → 1 we have been able to determine the complete cusp anomalous dimension A 4 for quarks at four loops. Our results for the n 0 f and n 1 f parts beyond the large-n c contribution are numerical, and lead to a relative accuracy of 10 −4 for these coefficients in QCD, which should be amply sufficient for phenomenological appli- article and of ref. [32] are in agreement with the form-factor calculations in refs. [46,47] for the planar part and refs. [117,118] for the n 2 f contributions. In particular, refs. [47] provided a quick confirmation of our result for the hardest part, the n f -independent contribution, and hence of our (mathematically not completely rigorous) reconstruction of the all-N expression of the large-n c non-singlet anomalous dimension.
The terms B 4 proportional to δ(1−x) in the four-loop splitting functions yield the universal eikonal anomalous dimension if properly combined with information on infrared singularities from the QCD form factor. This is an important ingredient for extending the threshold resummation for inclusive cross section to N 4 LL accuracy, i.e. (next-to-) 4 -leading logarithmic order.
In order to practically complete the QCD evolution equations at N 3 LO, corresponding results are required for the singlet splitting functions at four loops, i.e., the pure-singlet quark-quark splitting function P (3) ps (x) and those involving gluons, P qg (x), P gq (x) and P (3) gg (x). All these quantities are currently unknown beyond the N = 2 and N = 4 moments presented in ref. [31]. However, by following the approach of the present paper, it should be feasible to compute enough moments of the (theoretically much more complicated, see refs. [12,119]) corresponding flavour-singlet OMEs up to four loops with the Forcer program to gather sufficient information for first phenomenologically relevant approximations. We leave this topic to future research.
Incremental improvements in the flavour non-singlet sector can be obtained by calculating more moments, which is a hard problem within the present computational set-up for almost all colour factors, and by incorporating more external information, such as, e.g., a future exact result for the four-loop cusp anomalous dimension from calculations of the photon-quark form factor.
A derivation of the exact expressions for the n f -independent hardest parts will require, in addition, a much improved theoretical understanding. In this context it may be interesting to note that the ζ 5 part of γ (3)± ns (N ), which can be determined at all N from the presently available information (see appendix D) includes a contribution that vanishes in the large-n c limit. The resulting ln 2 N large-N behaviour needs to be compensated by non-ζ 5 terms, and it is tempting to identify the 5 ζ 5 in eq. (7.1) as the ζ 5 -'tail' of the function This function first occurred multiplied with positive powers of N in three-loop coefficient functions of DIS in ref. [22] and resurfaced, now multiplied with [S 1 (N )] 2 as in eq. (7.1), as the 'wrapping correction' in the anomalous dimensions in N = 4 maximally supersymmetric Yang-Mills theory [120], where it is crucial for obtaining the correct small-x limit, see, e.g., ref. [121]. We thus hypothesize that eq. (7.1) represents the first glimpse of the wrapping corrections in an anomalous dimension in QCD.
Form and Fortran files with our results can be obtained from the preprint server http://arXiv.org by downloading the source of this article. They are also available from the authors upon request. [115] to our notation. This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) grant MO 1801/1-2 and SFB 676 project A3, the European Research Council (ERC) Advanced Grant 320651, HEPGAME and the U.K. Science & Technology Facilities Council (STFC) grant ST/L000431/1. We also are grateful for the opportunity to use most of the ulgqcd computer cluster in Liverpool which was funded by the STFC grant ST/H008837/1. The Feynman diagrams have been drawn with the packages Axodraw [122] and Jaxodraw [123].

A Feynman rules
Below we present the Feynman rules for vertices arising from insertions of the operator O ns {µ 1 ,...,µ N } in eq. (2.1). All momenta p i are flowing into the operator vertex and we use where q is the outgoing momentum flow through the operator. The free Lorentz indices of the operator O ns {µ 1 ,...,µ N } are contracted with where the vector ∆ fulfils ∆ 2 = 0. We limit the derivation up to four additional gluons coupling to the operator, i.e., n = 6 in figure 1. For Feynman rules with up to three additional gluons and zero momentum flow through the operator, see also ref. [72] and references therein. The expressions for unpolarized quark operators in eqs. (A.2)-(A.6) are readily generalized to the polarized case by substituting / ∆ → / ∆γ 5 .

B Mellin moments at four loops
Here we present the anomalous dimensions at four loops for 1 ≤ N ≤ 16. Obviously, γ − ns (N = 1) = 0 at all orders. To fix our normalization, we write down the complete expression for N = 2 including all lower orders. Recall that a s = α s /(4π).

JHEP10(2017)041
The four-loop coefficient for the even moments N = 4 to N = 16 are given by

C Time-like splitting function
Here we present the difference between the space-and time-like non-singlet splitting functions at four loops, defined by δ P (3)± (x) ≡ P   A difference between the time-like and space-like case appears for the quantities P Here we finally present the exact expressions for the part of γ (3)± ns (N ) which is proportional to ζ 5 : Since only sums with w ≤ 2 can occur with ζ 5 , the corresponding functional form is so restricted that a direct determination and verification is possible with eight even and odd moments.