Large-nf Contributions to the Four-Loop Splitting Functions in QCD

We have computed the fourth-order nf^2 contributions to all three non-singlet quark-quark splitting functions and their four nf^3 flavour-singlet counterparts for the evolution of the parton distributions of hadrons in perturbative QCD with nf effectively massless quark flavours. The analytic form of these functions is presented in both Mellin N-space and momentum-fraction x-space; the large-x and small-x limits are discussed. Our results agree with all available predictions derived from lower-order information. The large-x limit of the quark-quark cases provides the complete nf^2 part of the four-loop cusp anomalous dimension which agrees with two recent partial computations.


Introduction
In the past years the next-to-next-to-leading order (NNLO) corrections in perturbative QCD have been determined for many high-energy processes, see Refs. [1][2][3][4][5][6][7][8] for some recent calculations. For processes with initial-state protons, NNLO analyses require parton distributions evolved with the three-loop splitting functions [9,10]. In some cases also the next-to-next-to-next-to-leading order (N 3 LO) corrections are important, e.g., for quantities with a slow convergence of the perturbation series or for cases where a very high accuracy is required. An example of the former is Higgs production at proton-proton colliders [11,12]. An example of the latter is the determination of the strong coupling constant α s from the structure functions F 2 and F 3 in lepton-nucleon deep-inelastic scattering (DIS), see Ref. [13], for which the N 3 LO coefficient functions have been obtained in Refs. [14,15]. In principle N 3 LO analyses of these processes require the four-loop splitting functions, although estimates of these functions via, for example, Padé approximants can be sufficient in some cases such as for DIS at large Bjorken-x.
At present a direct computation of the four-loop splitting functions P (3) ik (x) appears to be too difficult. Work on low-integer Mellin moments of these functions started ten years ago [16]; until recently only the N = 2 and N = 4 moments had been obtained of the quark+antiquark non-singlet splitting function P (3)+ ns together with the N = 3 result for its quark-antiquark counterpart P (3)− ns [17][18][19]. Using FORCER [20,21], a four-loop generalization of the well-known MINCER program [22,23] for the parametric reduction of self-energy integrals, it is now possible to derive more moments in the same manner as in Refs. [24][25][26] at the third order in α s . So far the moments up to N = 6 and N = 4 have been computed, respectively, for the non-singlet and singlet cases [27,28], and computations up to N = 8 are feasible. Further conceptual and/or computational developments are required, however, in order to obtain sufficient information for the construction of approximate x-space expressions analogous to those at three loops in Ref. [29].
The situation is far more favourable for the contributions to the functions P (3) ik (x) which are leading (in the singlet case) or leading and sub-leading (in the non-singlet case) in the number n f of effectively massless quark flavours. Here the harder four-loop diagram topologies do not contribute, and FORCER calculations above N = 20, and in some cases above N = 40, are possible. If suitably combined with information and expectations on the structure of these contributions in terms of harmonic sums [30,31], these fixed-N results turn out to be sufficient to find and validate the analytic dependence of these parts of the four-loop splitting functions on N, and hence on x in terms of harmonic polylogarithms [32], by LLL-based techniques [33][34][35]. This approach has been used before, e.g. in Refs. [36,37] for the three-loop transversity and helicity-difference splitting functions, and may be applicable to other four-loop quantities in the future. The present results include the n 2 f part of the four-loop cusp anomalous dimension also obtained in Refs. [38][39][40]. The remainder of this article is organized as follows: in Section 2 we set up our notations and briefly discuss the diagram calculations and the LLL analyses of the resulting integer-N moments. The analytic results for the n 3 f parts of P (3) ik and the n 2 f parts of P (3) ns in N-and x-space are presented and discussed in Sections 3 and 4. We summarize our results in Section 5.

Notations and calculations
The renormalization-group evolution equations for the dependence of the parton momentum distributions f a = u,ū, d,d, . . . , g of hadrons on the mass factorization scale µ f , form a system of 2n f +1 coupled integro-differential equations. These equations can be turned into ordinary differential equations by a Mellin transformation, and decomposed into 2n f −1 scalar (non-singlet) equations for the combinations of quark distributions and the 2×2 flavour-singlet quark-gluon system by using the general properties of QCD such as P gq i = P gq i = P gq . Note that P qg = 2n f P q i g .
The splitting functions in Eqs. (2.1) admit an expansion in powers of α s which we write as i.e., we identify (without loss of information) the mass-factorization and the coupling-constant renormalization scales. The difference between the splitting functions P + ns and P − ns for the first two non-singlet combinations in Eq. (2.4) and the pure-singlet quark-quark splitting function P ps = P qq − P + ns (2.6) starts at the second order in α s , the remaining difference at the third order in α s . To order α 4 s the latter quantity is proportional to the cubic group invariant d abc d abc /n c , while the other splitting functions can be expressed in terms of C F = 4/3 and C A = n c = 3 in QCD and quartic group invariants; the latter do not occur with the powers of n f that are considered in this article. The even-N or odd-N moments of the splitting functions are related to the anomalous dimensions γ(N) of twist-2 spin-N operators in the light-cone operator product expansion (OPE), see, e.g., Refs. [41,42]; we use the standard convention γ (n) (N) = −P (n) (N).
Our calculation of the four-loop splitting functions proceeds along the lines of Refs. [24][25][26]. The partonic DIS structure functions are mapped by the optical theorem to forward amplitudes probe (q) + parton(p) −→ probe (q) + parton (p) (2.8) with p 2 = 0 and q 2 = −Q 2 < 0. Via a dispersion relation their coefficients of (2p · q/Q 2 ) N then provide, depending on the structure function under consideration, the even-N or odd-N moments of the unfactorized partonic structure functions. These quantities are calculated in dimensional regularization with D = 4 − 2ε, and the n-loop splitting functions can be extracted from the coefficients of ε −1 α n s . For the even-N determination of the splitting functions P + ns and P ik in Eq. (2.4) we use the photon and the Higgs boson in the heavy-top limit as the probes. The splitting functions P − ns and P v ns are determined from the odd-N vector -axial-vector interference structure function F 3 . The projection on the N th power in the parton momentum p leads to self-energy integrals that can be solved by the FORCER program. The complexity of these integrals increases by four if N is increased by two. Together with the steep increase of the number of integrals with N, see the discussion of the harmonic projection in Ref. [23], this limits the number of moments that can be calculated. So far high values of N cannot be reached for the top-level 4-loop diagram topologies.
The raw diagram databases provided by QGRAF [43] are heavily manipulated by (T)FORM [44][45][46] programs to provide the best possible starting point for the main integral computations. As discussed in Ref. [47], one important step is the identification of ℓ-loop self-energy insertions, which reduces many n-loop diagrams to fewer (n − ℓ)-loop diagrams in which one or more propagators have a non-integer power. For the large-n f contributions under consideration in the article, genuine four-loop diagrams remain after this step only in the calculation of the C A n 3 f part of P (3) qg , and these diagrams have a rather simple topology: in the notation of MINCER they are generalizations of the Y3 and O1 three-loop topologies. The hardest diagrams occur in the C A C F n 2 f and n 2 f d abc d abc /n c non-singlet cases: these are three-loop BE topologies with a one-loop gluon propagator, see Fig. 2; the highest N calculated here for any of these is N = 27.
As far as they are known from fixed-order calculations [9,10] and all-order resummations of leading large-n f terms [48][49][50], the even-N or odd-N moments of the splitting functions (i.e. the anomalous dimensions) can be expressed in terms of simple denominators, D k a = (N + a) −k and harmonic sums [30,31] with 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 = 2n − 1 occur in the n-loop anomalous dimensions, but no sums with an index −1.
For terms with D k a and/or coefficients that include values ζ m of the Riemann ζ-function (with m ≥ 3, ζ 2 does not occur in these functions), the maximal weight of the sums is reduced by k + m. It is, of course, possible that other structures occur in the n-loop anomalous dimensions at n ≥ 4 -already the three-loop DIS coefficient functions include terms where special combinations of sums are multiplied by low positive powers of N [14,15]. However, one may expect this to happen at n = 4 only in the terms with low powers of n f which receive contributions from generically new diagram topologies. Disregarding new structures and terms with ζ m≥3 which are much easier to fix from low-N results, a general ansatz for the n-loop anomalous dimensions then is where S w (N) is a shorthand for all harmonic sums with weight w and S 0 (N) ≡ 1. The terms with c 00w only occur in the quark-quark and gluon-gluon splitting functions and are restricted by the known large-N structure of these functions [51][52][53]. In all cases the range of the sums is reduced for large-n f contributions in a manner that can be inferred from the results at n ≤ 3 and from the prime-factor decompositions of the denominators of the calculated moments.
Even so, Eq. (2.11) usually includes far too many coefficients for a direct determination from as many calculated moments. These coefficients, however, are integer modulo some predictable powers of 1/3 at n ≤ 2 [9,10] and in Refs. [48][49][50]. Hence the systems of equations can by turned into Diophantine systems which require far fewer equations than unknowns. Given the present limitations of the calculation of diagrams with BE topology, this is still not sufficient for the n 2 f contributions to the four-loop non-singlet splitting functions. However, these functions include additional structures that facilitate solving these equations with the calculable moments.
The crucial point for the determination of the n 2 f parts of γ (3)± ns (N), already presented in [27], is to write its colour-factor decomposition in two ways, where the result has been rendered more compact by using the abbreviations As the overall leading-order quantity P (0) qq , the splitting function corresponding to Eq. (2.14) is the same for the present initial-state and the final-state (fragmentation distributions) evolution, cf. Refs. [54][55][56], and invariant under the x-space transformation f (x) → x f (1/x). The (combinations of) harmonic sums in Eq. (2.14) are 'reciprocity respecting' (RR), i.e., their Mellin inverses are invariant under the above transformation. The same holds for the combinations of denominators in Eq. (2.15). Except for S 2 1 and S 3 1 -products of RR sums lead to higher weight RR sumsall reciprocity-respecting sums to weight three occur in Eq. (2.14). The list of RR function to this weight has been given in Ref. [36] with a slightly different basis choice at w = 3.

Like the overall NLO anomalous dimensions γ
(1)± ns (N), the next-to-leading order d abc d abc /n c contribution γ (3)s ns (N) is not reciprocity-respecting. However, and this is the crucial point, its RRbreaking part can be calculated from Eq. (2.14) according to the conjecture of Ref. [53]. For the n 2 f contribution addressed here it is given by − 2 3 n f d dN γ (2)s ns (N), where the differentiation can be carried out, for example, via the asymptotic expansion of the sums, see also Ref. [57]. That leaves an unknown reciprocity-respecting generalization of the form (2.14) with additional w = 4 sums which can be chosen as Including also ν 2 terms, one arrives at a trial function with 79 coefficients, of which as many as 15 can be eliminated by imposing the existence of the first moment and the correct values (zero) for its ζ-function contributions, and 9 can be assumed to vanish (all contributions with S 3 1 and S 4 1 ). The remaining 56 coefficients have then been found using the 12 odd moments with 3 ≤ N ≤ 25.
The correctness of the solution has been verified by the (non-ζ) value of the first moment and the result at N = 27. It is possible, though, to judge 'by inspection' whether a solution returned by the Diophantine equation solver [34,35] is correct. For example, the above solution is returned as A pattern such as the one above for the about 30 coefficients of the highest-weight functions, with larger and more random coefficients at the left (low-weight) end, is a hallmark of a correct solution. In fact, correct and incorrect solutions were correctly identified by inspection in all present calculations as well as in the preparation of Ref. [37].
Of the n 3 f contributions to the singlet splitting functions in Eq. (2.4), only the case of P qg is critical. Unlike the other three cases this function is suppressed by only two powers of n f relative to the lowest-n f term, recall the remark below Eq. (2.4), and includes contributions from sums up to weight four instead of weight three. Hence a considerably larger basis set is required in Eq. (2.11). At the same time the fixed-N calculations are harder for P (3) qg than for the other three cases, in particular for the C A n 3 f contribution, as already indicated on p. 4. Yet, using reasonable assumptions based on the three-loop splitting function, we managed to find suitable functional forms with 101 unknown coefficients for the C F n 3 f part (with only positiveindex sums but overall weight up to six) and 115 unknown coefficients for the C A n 3 f part (including alternating sums but an overall weight of five), which we were able to determine from the even moments 2 ≤ N ≤ 40 in the former and 2 ≤ N ≤ 44 in the latter case. Several higher moments were employed for the validation of the C F n 3 f result and the C A n 3 f coefficients were checked using N = 46. Some of the four-loop and three-loop C A n 3 f diagrams at N > 40 were calculated using an alternative approach for generalized Y and O MINCER topologies that avoids the harmonic projection [23]. This approach may be reported on later in a more general context.

Results in N-space
In this section we present the analytic expressions for the n 2 f and n 3 f contributions to the three nonsinglet anomalous dimensions and the n 3 f parts of their four flavour-singlet counterparts in the MS scheme. As in Eqs. (2.13) and (2.14) above, all harmonic sums (2.9) and (2.10) have the argument N which is suppressed in the formulae for brevity.
The results for γ (3)± ns are presented in terms of the decomposition (2.12). The large-n c part is the same for these two cases, while the contributions with the 1/n c -suppressed 'non-planar' colour factor (C A − 2C F ) are valid at even N for B As for the complete corresponding three-loop quantities in Ref. [10], the difference between the odd-N result (3.3) and the even-N result (3.2) is much simpler than those expressions and given by (3.5) The leading large-n f contribution is the same for the three types of non-singlet quark distributions in Eq. (2.3). It has been obtained to all orders in α s in Ref. [48]. Our results agree with the corresponding fourth-order coefficient which in our notation reads, for even and odd N, The new functions (3.1) -(3.5) are illustrated in Fig. 2. The results at non-integer values of N have been calculated by a numerical Mellin transformation of the x-space expressions in the next section; for the analytic continuation in N of the harmonic sums to weight five see also Ref. [57].
Up to terms suppressed by two powers of 1/N, also the large-N behaviour of the three nonsinglet anomalous dimension is the same with where γ e is the Euler-Mascheroni constant. The coefficients A n are relevant beyond the evolution of the parton distributions, since they are identical to the n-loop cusp anomalous dimensions [51]. The result at three loops can be found in Eq. (3.11) of Ref. [9], its n f part was derived before in Refs. [58,59]. Our new results (3.1) and (3.2) specify the n 2 f coefficient of A 4 . Together with the long-known n 3 f result [48,60] given by the large-N limit of Eq. (3.6) we obtain The large-n c limit of this result has also been derived in Ref. [38], and the C 2 F n 2 f part in Refs. [39,40]. Hence all n 2 f contributions in Eq. We now turn to the leading large-n f anomalous dimensions for the even-N flavour-singlet evolution (2.4), starting with the pure singlet contribution (2.6): (3.10) As expected from the lower orders, the highest-weight sums in the four-loop off-diagonal contributions are proportional to the leading-order structures Using these abbreviations, the fourth-order leading-n f parts of the gluon-quark and quark-gluon anomalous dimensions are given by  Finally the corresponding contribution to the gluon-gluon anomalous dimension reads (3.14) The C A part of Eq. (3.14), which is a non-singlet -type quantity and hence could be written in a more compact manner in terms of the quantities in Eq. (2.15), has been obtained already in Ref. [50]. Its leading large-N coefficient is related to that in Eq. (3.8) by the 'Casimir scaling' C A /C F . Moreover two linear combinations of Eq. (3.13) with Eq. (3.10) and the C F part of Eq. (3.14) were derived in Ref. [49,50]; our results agree also with those findings. Eq.

Results in x-space
The four-loop splitting functions P (3) ik (x) are obtained from the above N-space results by an inverse Mellin transformation which expresses these functions in terms of harmonic polylogarithms. This transformation can be performed by a completely algebraic procedure [32,62] based on the fact that harmonic sums occur as coefficients of the Taylor expansion of harmonic polylogarithms.
Before we present our results, we recall the basic definitions [32]: The lowest-weight (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 where we have used the abbreviation order results note the discussion below Eq. (30) in Ref. [63]. Finally the common leading large-n f contribution to the N 3 LO evolution of all three types of quark distributions in Eq. (2.3) reads P (3) ns (x) Also the large-x limit is the same for the three non-singlet splitting functions. It is given by in terms of the same constants as in Eq. (3.7), i.e., the n a >1 f contributions to A 4 and C 4 have been given in Eqs. (3.8) and (3.9). The coefficients B 4 can be read of from Eqs. (4.6), (4.7) and (4.12). The difference between P − ns and P + ns and the splitting function (2.7) are suppressed by two powers of (1−x) with respect to the leading term in Eq. (4.13).
The non-singlet splitting functions include double-logarithmic small-x contributions up to ln 2ℓ x at N ℓ LO. The coefficients of these leading-logarithmic (LL) parts of P ± ns have long been known to all orders [64,65]; the presence of a ln 4 x term in P s ns at NNLO had not been predicted before the three-loop calculation in Ref. [9]. Contributions where k powers of (C A , C F ) in the colour factor are replaced by n k f are suppressed by k powers of ln x relative to the overall leading logarithms. Hence we expect terms up to ln 4 x and ln 5 x, respectively, in P  (4.17) The analytic structure of the LL resummations is very different for P + ns and P − ns with [64,65] and where z = N(2a s N c ) −1/2 , and D p (z) denotes a parabolic cylinder function [66]. The expansion of Eq. (4.19) in powers of a s is an asymptotic expansion, in contrast to Eq. (4.18). The difference between the two expansions vanishes in the large-n c limit.
An extension of these resummations to next-to-leading logarithmic (NLL) accuracy and beyond is known so far only for the former case -for the x 2n ln ℓ x terms at n ≥ 0; for the x 2n+1 ln ℓ x terms the roles of P + ns and P − ns are interchanged in this respect. A determination of the N ℓ LL terms on the basis of N ℓ LO information is possible from the D-dimensional structure of the unfactorized expressions, analogous to the case of the final-state splitting functions and coefficient functions in semi-inclusive annihilation [67,68]. The first term in Eq. (4.15), an overall NNLL contribution, agrees with the result in Eq. (4.6) of Ref. [69] after α s -expansion and Mellin inversion. For details and results on the singlet cases and coefficient functions see Ref. [70].
The n 2 f contributions to the three non-singlet splitting functions are illustrated In Figs. 5 and 6 on linear and logarithmic scales in x. The latter have been extended to x = 10 −6 in order to include the onset of the steep small-x rise of all these functions. The difference (4.10) between the n 2 f parts of P At asymptotically small values of x, the behaviour of these functions is given by their respective leading ln 4 x and ln 5 x logarithms in Eqs. (4.14) -(4.16). As shown in the figures, though, the onset of the resulting steep rise towards x = 0 is delayed to x ≈ 10 −5 by the effect of the non-leading logarithms. In fact, even at the lowest x-values shown here a relevant approximation for P + ns and P − ns is obtained only if all ln x terms are taken into account. The situation is more favourable for P s ns but, unlike for the three-loop contribution [9] to this function, also here the leading logarithmic result is totally different from the actual function at all physically sensible values of x.  Figure 6: As the right panel of Figure 5, but for the splitting functions P − ns (left) and P s ns (right). Due to Eq. (2.5) all numbers have to be divided by (4π) 4 ≃ 25000 for an expansion in powers of α s . The x-space splitting functions corresponding to the flavour-singlet anomalous dimensions in Eqs. (3.10) -(3.14) are given by Like their non-singlet counterparts, the singlet splitting functions receive a double-logarithmic small-x enhancement of the form α n s ln ℓ x with 0 ≤ ℓ ≤ 2n. However, the small-x behaviour in the singlet case is dominated by additional single-logarithmic x −1 ln ℓ x terms, see Refs. [77][78][79][80]. In the present α 4 s n 3 f cases, only non-logarithmic x −1 terms occur and the small-x expansions read The coefficients of ln 4 x in these results agree with the results of the double-logarithmic small-x resummation [70]. The pattern in Eq. (4.30), no small-x logarithms, is the same as for the C F n 2 f contribution to P gq at order α 3 s . The x −1 terms in Eqs. (4.29) and (4.31) show an interesting feature in the large-n c limit C F → 1 2 n c : the resulting coefficients of x −1 n c for P (3) qg and P (3) gg are identical to those of x −1 C F for P (3) ps and P (3) gq in Eqs. (4.28) and (4.30), respectively. For the QCD values of the colour factors, the ratio between the x −1 coefficients is 2.11 for the upper-row splitting functions and 2.09 for their lower-row counterparts; hence these ratios are between their overall large-n c limit of 2 and the Casimir-scaling value of 9/4.

Summary
As a first step towards the determination of the N 3 LO splitting functions P (3) ab (x) in perturbative QCD beyond the leading large-n f results of Refs. [48][49][50], we have derived the complete n 2 f parts of the four-loop non-singlet quark-quark splitting functions and all n 3 f contributions to their flavoursinglet counterparts in the MS scheme. These results have been obtained by analytically computing a fairly large number of Mellin moments N in the approach of Refs. [24][25][26] -made possible by the development of the FORCER program [20,21] for the computation of massless four-loop selfenergy integrals -and a subsequent determination of the all-N and all-x expressions using the number-theoretical results and tools of Refs. [33][34][35], a method that has been applied already to three-loop splitting functions in Refs. [36,37].
Our results agree with Refs. [48][49][50], with the pioneering low-N non-singlet computations of Refs. [17][18][19], and with the recent determinations of n f contributions to the four-loop cusp anomalous dimension [38][39][40] which appear in our results as the coefficient of ln N at large N or 1/(1 − x) + in the large-x expansion. We also agree with the prediction of Ref. [53] for the coefficient of ln (1−x) in the non-singlet cases and, in the small region of overlap, with the resummations of highest three small-x and large-x double logarithms in Refs. [69,70,74,75]. Most interestingly our results are in agreement with the remarkably simple (if incomplete -the ζ 2 (C A − 2C F ) contributions are excluded) generalization of the leading-log small-x resummation [64,65] for the quark+antiquark non-singlet splitting function P + ns to all powers of ln x proposed in Ref. [71]. By themselves the present results are not phenomenologically useful. We hope, though, that it will be possible to complement them in the near future by approximate expressions of the remaining (and numerically more important) contributions to the functions P (3) ab (x), analogous to those employed at NNLO [29] before the results [9,10] became available, and hence facilitate improved N 3 LO analyses of DIS and hard processes at colliders. One may also hope that the present results will provide useful additional 'data' for future studies of the structure of the perturbation series for the splitting functions which, in turn, may lead to more explicit four-loop calculations and results.
FORM [44][45][46] files of our N-space expressions in terms of harmonic sums [30,31] and their x-space counterparts in terms of harmonic polylogarithms [32] can be obtained from the preprint server http://arXiv.org by downloading the source of this article. Furthermore they are available from the authors upon request.