Resummation of small-x double logarithms in QCD: inclusive deep-inelastic scattering

We present a comprehensive study of high-energy double logarithms in inclusive DIS. They appear parametrically as alpha_s^n ln^{2n-k} x at the n-th order in perturbation theory in the splitting functions for the parton evolution and the coefficient functions for the hard scattering process, and represent the leading corrections at small $x$ in the flavour non-singlet case. We perform their resummation, in terms of modified Bessel functions, to all orders in full QCD up to NNLL accuracy, and partly to N^3LL and beyond in the large-n_c limit, and provide fixed-order expansions up to five loops. In the flavour-singlet sector, where these double logarithms are sub-dominant at small x compared to single-logarithmic alpha_s^n x^{-1} ln^{n-k} x BFKL contributions, we construct fixed-order expansions up to five loops at NNLL accuracy in full QCD. The results elucidate the analytic small-x structure underlying inclusive DIS results in fixed-order perturbation theory and provide important information for present and future numerical and analytic calculations of these quantities.


Introduction
Inclusive deep-inelastic lepton-hadron scattering (DIS) is an experimental and theoretical reference process for Quantum Chromodynamics (QCD), the theory of the strong interaction. Important information on the parton (quark and gluon) distribution functions (PDFs) of the proton, in particular, is provided by the dependence of the corresponding cross sections or the structure functions F a (x, Q 2 ) on the Bjorken variable x and on the scale Q 2 = −q 2 set by the momentum q of the exchanged (gauge) boson. Moreover the scaling violations, i.e., the Q 2 -dependence of the structure functions F 2 and F 3 , facilitate high-precision determinations of the strong coupling constant α s .
Due to their relation to propagator-type Feynman integrals via forward Compton amplitudes and the light-cone operator-product expansion (OPE), see, e.g., refs. [1][2][3][4] and references therein, structure functions in DIS are particularly well suited for analytical high-order computations in massless perturbative QCD. Indeed, the complete third-order contributions to the (initial-state) splitting functions governing the evolution of the PDFs were obtained more than fifteen years ago [5,6] in computations that also provided the third-order cross-section projections (coefficient functions) for the most important structure functions in spin-averaged DIS [7][8][9]. During the past five years those computations have been extended, if only for a limited number of Mellin moments, to the fourth order in α s [10][11][12][13]; for the lowest moments see also refs. [14][15][16][17][18].
The perturbation series for the splitting functions and coefficient functions appear to be very well behaved (except for the longitudinal structure function F L [7,10]) outside the threshold region 1−x ≪ 1 and the high-energy limit x → 0. With the exception of the diagonal (quarkquark and gluon-gluon) splitting functions in the standard MS scheme [19][20][21], the splitting and coefficient functions include threshold double logarithms at all powers of 1−x. The dominant (1−x) −1 ln ℓ (1−x) contributions to the coefficient functions in DIS have been resummed to a high accuracy [22,23], see also ref. [24], in the framework of the soft-gluon exponentiation [25][26][27][28][29][30]. The resummation of the double logarithms has been extended to non-negative powers of 1−x by analyzing the physical evolution kernels of the structure functions [31][32][33], see also ref. [34], and the structure of the 'raw' (unfactorized) expressions in dimensional regularization [35][36][37].
The latter approach can be applied also to high-energy double logarithms ∼ α n s x p ln 2n−n 0 −k x in splitting functions and coefficient functions, albeit, at least in its present form, not at all powers p of x. In particular, the resummation of the dominant p = −1 contributions to the splitting functions for the final-state parton fragmentation functions and to the coefficient functions for semiinclusive electron-positron annihilation (SIA) [38,39] have been extended to the k = 2 next-tonext-to-leading logarithmic (N 2 LL) accuracy in refs. [40,41], see also refs. [42,43]. Corresponding p = 0 results for inclusive DIS were obtained at about the same time. While being formally analogous to their SIA counterparts, these results were not of direct phenomenological relevance, and only one example expression was presented at the time [44].
Such results become relevant, however, for approximate or exact reconstructions of higherorder splitting functions and coefficient functions, if they can be combined with a sufficient amount of other information, such as a sufficiently large number of Mellin moments. Due to the develop-ment of the FORCER program [45] for four-loop propagator-type Feynman integrals, this point has now been reached for fourth-order corrections; see refs. [11,46,47] for published results on the splitting functions. In fact, in the latter two articles fourth-order predictions of the resummations discussed above and of a complementary proposal of ref. [48] have already been employed and, where feasible, confirmed. Hence it is now timely to present, in sufficient detail to assist future research, the status of the resummation of small-x double logarithms for inclusive DIS.
The remainder of this article is organized as follows: In section 2 we specify our notation and discuss the available formalisms, and their limitations, to the resummation of small-x (double) logarithms. We also briefly indicate how the calculations have been performed. The results for the splitting function for the evolution of flavour differences of sums of quark and antiquark PDFs are presented in section 3. This is the case for which the two complementary approaches overlap. The results include another striking illustration of the phenomenological inadequacy of representing the splitting functions, at any relevant x, solely by a N ℓ LL small-x approximation at some fixed ℓ.
In sections 4 and 5 we present the N 2 LL predictions for the corresponding non-singlet coefficient functions and for the flavour-singlet splitting and coefficient functions. In view of the findings in section 3, we focus in these sections on fourth-and fifth-order predictions and the allorder structure of the x 0 ln k x contributions. We expect that the former results will become useful in combination with large-x information on these functions, while the leading-logarithmic all-order expressions may provide useful 'data' for future research into the small-x structure of splitting functions and coefficient functions in DIS. We briefly summarize our findings in section 6. Some additional material that may be useful to future research can be found in the appendix.

Notation, formalism and calculations
Disregarding 1/Q 2 power corrections, the structure functions in DIS can be generically written as F a (x, Q 2 ) = C a,i (a s ) ⊗ f i (Q 2 ) (x) (2.1) in terms of the coefficient functions C a,i (x, a s ) and the corresponding (combinations of) parton distributions f i (x, Q 2 ). Here and below we identify the renormalization and mass-factorization scales µ 2 R and µ 2 F with the physical scale Q 2 ; the dependence on µ 2 R and µ 2 F can be readily reconstructed a posteriori, see, e.g., sections 2 of refs. [49,50]. ⊗ represents the Mellin convolution, given by The coefficient functions C a in eq. (2.1) and the splitting functions P ik in eq. (2.3) can be expanded in powers of the strong coupling constant, which we normalize as a s = α s (Q 2 )/(4π), P (x, a s ) = a s P (0) (x) + a 2 s P (1) (x) + a 3 s P (2) (x) + a 4 s P (3) (x) + . . . , (2.4) C a (x, a s ) = c (0) a (x) + a s c (1) a (x) + a 2 s c (2) a (x) + a 3 s c where i = q, g. Consequently, the terms up to c (n) and P (n) form the (next-to) n -leading order (N n LO) approximation of perturbative QCD for F 2 and F 3 , while c (n+1) L and P (n) are required for this accuracy for the longitudinal structure function F L .
The coefficients in eqs. (2.4) and (2.5) include high-energy double logarithms, with contributions up to x p ln 2n x for P (n) and c (n+1) L , and terms up to x p ln 2n+1 x for c (n+1) 2 and c (n+1) 3 at p ≥ 0 1 . Our main approach to the resummation of these logarithms, i.e., to the determination of the coefficients of a n s x p ln 2n−n 0 −k x contributions to eqs. (2.4) and (2.5), is analogous to that presented in ref. [40] for the case of p = −1 in semi-inclusive e + e − annihilation.
The primary objects of this resummation are the unfactorized partonic structure functions in dimensional regularization (we use D = 4 − 2ε) where, for simplicity, the indices labelling different structure functions and parton distributions have been suppressed. The functions C are given by Taylor series in ε; the ε m terms c (n,m) a including m more powers in ln x than the 4-dimensional (ε = 0) coefficient functions (2.5). The transition functions Z consist only of negative powers of ε and can be written in terms of the splitting functions (2.4) and the expansion coefficients β n of the D-dimensional beta function, with [54][55][56][57] etc where C A = n c = 3 and C F = (n 2 c − 1)/(2n c ) = 4/3 in QCD. Here and below, n f denotes the number of light flavours.
In Mellin N-space, the transition functions are related to the splitting functions by (2.10) Using this relation, Z(N, a s ) can be expressed order-by-order in terms of the anomalous dimensions γ n (N) = −P (n) (N) and β n . The first four orders in a s , allowing for Z and γ n to be matrices, read The higher-order contributions have been generated using FORM and TFORM [58][59][60] to a sufficiently high order for the computations of this paper. At order a n s , the dependence of Z on β m and γ m can be summarized as ε −n : γ 0 , β 0 , ε −n+1 : . . . , γ 1 , β 1 , ε −n+2 : . . . , γ 2 , β 2 , . . . , ε −1 : γ n−1 . (2.12) Hence fixed-order knowledge at N m LO (i.e., of the splitting functions to P (m) , beta function to β m and the corresponding coefficient functions) fixes the first m+1 coefficients in the ε-expansion of F at all orders in a s . Furthermore, the property P (n) ∼ ln 2n x ⇔ γ n ∼ N −2n−1 means that β 0 and β 2 0 enter eq. (2.11) at the next-to-leading logarithmic (NLL) and N 2 LL level, while β 1 contributes only from the fourth (N 3 LL) logarithms. Similarly, β 2 enters only at N 5 LL accuracy and beyond.
If the above N m LO knowledge can be extended to all powers in ε at a given logarithmic accuracy, e.g., to all coefficients of a n s ε −n+k N −n−k+ℓ or a n s ε −n+k ln n+k−1−ℓ x for F 2 in eqs. (2.7) and (2.11) for a fixed ℓ, then we arrive at an all-order resummation of these terms, e.g., of the N ℓ LL x 0 contributions to the splitting functions and coefficient functions contributing to F 2 .
This situation is completely analogous to that of the large-x double logarithms in refs. [35][36][37]. In that case, the structure that allows the extensions to all ε has been inferred from the calculations of inclusive DIS via suitably projected gauge-boson parton cross sections as carried out at two loops in refs. [61][62][63]. The same strategy can be applied here. The maximal (2 → n + 1 particles) phase space for these processes at order α n s can be schematically written as [64,65] 1 If the integrals for the n-th order purely real (tree graph) contributions F do not lead to any further factors x ε , their expansion (for a = L) around x = 0 can be written as If, furthermore, the mixed real-virtual contributions (2 → r + 1 particles with n − r ≥ 1 loops) include no more than n − r additional factors of x ε from the loop integrals, then we arrive at a,i,p + . . . (2.15) with n a = 1 and l a = n for all structure functions considered here, with the exception of F L for which n a = 3 and l a = n − 1. By expanding x εl in powers of εl ln x, it can be seen that A, B and C are the coefficients of the LL, NLL and NNLL (≡ N 2 LL) contributions, respectively. For a given value of p, the N-space counterpart of eq. (2.15) reads Since eq. (2.11) includes only poles up to ε −n at order a n s , the terms with ε −2n+1 , . . . , ε −n−1 have to cancel (for a = L) in the sums (2.15) and (2.16). Hence there are n−1 'zero' relations between the LL coefficients A (n,l) , n−2 relations between the NLL coefficients B (n,l) etc. Moreover, as mentioned above, the N m LO results provide the (non-vanishing) coefficients of ε −n , . . . , ε −n+m at all orders n, and thus m + 1 additional relations between the coefficients in eqs. (2.15) and (2.16). Consequently the highest m+1 double logarithms, i.e. the N m LL approximation, can be determined and, except for the N m LL terms at order a m+1 s , over-constrained order-by-order from the N m LO results. This feature also holds for a = L but (due to c (0) L,i (x) = 0) with only n−2 'zero' relations but also one term fewer in the above sums.
Using the known unfactorized N 2 LO expressions for the structure functions, it is now possible to establish, to all orders, for which cases eqs. (2.15) and (2.16) hold. It turns out that these equations, and hence the resulting resummation of small-x double logarithms, are applicable at Here 'e.m.' (electromagnetic) denotes photon exchange, and F φ is the structure function for DIS via the exchange of a scalar that, like the Higgs boson in the heavy-top limit, couples directly only to gluons. The third-order coefficient functions for this structure function, which is experimentally irrelevant but theoretically useful, have been presented in ref. [33]; the second-order results have also been obtained in ref. [66]. For the coefficient functions for ν −ν charged-current DIS see refs. [67][68][69]. For completeness, we have included the most important structure function g 1 in spin-dependent DIS, for its coefficient functions and splitting functions at NNLO see refs. [70][71][72][73][74].
It is worthwhile to note that all structure functions in eq. (2.17) are accessible only at even N via forward Compton amplitudes or the OPE; conversely all structure functions in eq. (2.18) are odd-N based (for a detailed discussion see, e.g., ref. [67]) -a fact that can hardly be a mere coincidence. Moreover, the differences between the splitting functions and coefficient functions for the flavour non-singlet ν +ν and ν −ν structure functions vanish in the limit of a large number of colours n c , see refs. [5,46,47,68,75].
The structural difference between these two cases, and its suppression at large n c , can already be seen from the p = 0 LL resummation of the corresponding even-N and odd-N splitting functions P + ns and P − ns for the flavour differences of quark-antiquark sums and differences [76][77][78]: where z = N (2 a s n c ) −1/2 , and D p (z) denotes a parabolic cylinder function [79]. Note that the expansion of eq. (2.20) in powers of a s is an asymptotic series, in contrast to eq. (2.19).
In ref. [48] a surprisingly simple generalization has been proposed of the equation, first derived in ref. [76], that leads to eq. (2.19). This generalization can be stated as up to terms that are large-n c suppressed and include even-n values ζ n of Riemann's ζ-function. 2 Inserting the Laurent expansion about N = 0 into eq. (2.21), one can readily solve this relation to 'any' desired order n for the coefficients p + n,k with k ≤ 2n − 1 that correspond to powers of ln x in the small-x expansion. Specifically, the coefficients p + n,0 to p + n,2m+1 can be predicted (with the above restriction) at all n > m from eq. (2.21) with k max = 2m + 1 if P + ns (N) is completely known to N m LO. So far these predictions have been verified for the n 2 f and n 3 f terms and the complete large-n c limit at N 3 LO [46,47].
The above two approaches to the resummation of small-x terms overlap for the x 0 ln k x double logarithms of P + ns (x, α s ), but are largely complementary otherwise. The predictions of eq. (2.16) cover far more than just x 0 part of P + ns (x, α s ), while eq. (2.21) is very powerful in this specific case, in particular in the large-n c limit.
Both approaches require Laurent expansions of the fixed-order input quantities, including nonnegative powers of N + p, as written down at p = 0 for P (n)+ ns (N) in eq. (2.22). These expansions can be obtained, for example, by expanding the exact x-space expressions in terms of harmonic polylogarithms (HPLs) [83], using the HARMPOL package for FORM [58] together with M x p ln k 1 (2.23) An easy extension to the coefficients of non-negative powers of N+p is to transform the functions to N-space harmonic sums [84,85], multiply by a sufficiently large power s of 1/(N+p), transform back to x-space, proceed as above, and finally multiply by (N+p) s . Routines for the Mellin transform of the HPLs and its inverse are also provided by the HARMPOL package. For the convenience of the reader, the p = 0 coefficients employed in this article are collected in appendix A.
The resummation predictions for the splitting and coefficient functions can then be computed order by order in α s . Using FORM and TFORM [58][59][60], this has been done up to order α 30 s and α 60 s , respectively, for the flavour singlet and non-singlet cases. Using the formal similarity to the SIA cases covered in refs. [40,41], these results can then be employed to infer their generating functions via over-constrained systems of linear equations, thus arriving at all-order expressions.

Results for the non-singlet splitting functions
Non-singlet quantities are dominated at small x by their x p ln k x contributions with p = 0 and k ≥ 0 which correspond to poles at N = 0 in Mellin space. These terms can be resummed via eq. (2.16) for the splitting function P (n)+ ns which enters the structure functions in eq. (2.17). The resummed N-space expressions can be expressed in terms of (3.1) The N 2 LL result for P (n)+ ns , already presented in ref. [44], can be written as where β 0 in eq. (2.8) has been used instead of n f for a more compact representation. The two terms in the first line of eq.
and P +(4) to the moments of eq. (2.4). The n 2 f part of eq. (3.3) and its complete large-n c limit have been employed in refs. [46,47], respectively, as constraints in the determination of the all-N expressions from a limited number of moments and endpoint constraints. Conversely, the verification of those all-N expressions -via results at higher N and independent form-factor calculations [86][87][88][89][90] that include the large-N limit of P ns , the (light-like) four-loop cusp anomalous dimension [19,20] provides a stringent check of eq. (3.3).
Using the N 3 LO results [46,47] together with the corresponding coefficient function for F 2 in ref. [8], eqs. (3.2) -(3.4) can be extended to N 3 LL small-x accuracy for the next-to-leading large-n f terms and in the large-n c limit. The latter results will be presented below, together with the predictions of eq. (2.21) at this and higher orders in the small-x expansion.
The structure of the closed-form N-space expression (3.2) is similar to the 'non-singlet' p = −1 part of the 'time-like' splitting function P T gg (N) for final-state fragmentation functions. The crucial difference is the sign of ξ in eq. (3.1) which leads to qualitative differences. In x-space, the latter splitting function can be expressed [41] in terms of Bessel functions which exhibit an oscillatory behaviour in the small-x limit. In fact, the resummation is found to completely remove the huge small-x spikes present in the fixed-order results for the time-like splitting functions [91,92].
In the present 'space-like' (initial-state) case, on the other hand, the N 2 LL x-space expression is given by and in terms of modified Bessel functions I n (z), see section 9.6 of ref. [93]. The first two terms in the curly bracket in the first line are the LL and NLL results, respectively; the remaining terms provide the N 2 LL contribution. The latter can be written in different ways due to the recurrence relation expressing I n+1 in terms of I n and I n−1 . The form chosen above yields the most compact coefficients and is in line with our basis choice for the higher-accuracy large-n c expressions below.
The functions I n (z) have the exponential form (2π z) −1/2 e z 1 + O(z −1 ) in the large-z limit, and thus dwarf the fixed-order small-x behaviour in a very unstable manner: the LL result is positive, the NLL contribution is negative (for the physically relevant case β 0 > 2C F ), the N 2 LL 'correction' is positive etc. It is interesting to note that the coefficients of the (for x → 0 at fixed a s ) asymptotically dominant a s (a s ln x) ℓ I 1 (z) terms in the first line of eq. (3.5), which correspond to the S −2ℓ+1 parts of the N ℓ LL contributions in eq. (3.2), seem to point towards a 'second resummation', a feature already observed in ref. [41] for the time-like splitting functions. We will return to this issue with more 'data' below, see eqs. (3.16) and (3.17).
As discussed in section 2, the contributions up to N 5 LL -with the exception of large-n c suppressed terms containing ζ 2 -can be obtained from eq. (2.21) [48] and P + ns up to three loops [5]. The resulting coefficients at N 3 LO (four loops), defined in eq. (2.22), are given by up to large-n c suppressed contributions with (powers of) ζ 2 which may include quartic group invariants. The large-n c limits of eqs. (3.8) -(3.10) have been employed (and verified, recall the discussion below eq. (3.4)) together with that of eq. (3.3) in ref. [47]. The n 2 f terms agree with eq. (4.14) of ref. [46] which, unlike eqs.  [94,95], agree with eq. (4.17) of ref. [46].

The corresponding N 4 LO (five-loop) coefficients read
ζ 2 n f n 3 c C F + 11776 27 The large-n c limit of eq. (3.11) can also been also obtained via eqs. (2.7) to (2.16) using the N 3 LO result P +(3) ns of ref. [47].
Using the all-N large-n c expression for P ns , it is possible to predict also the N −3 and N −2 (N 6 LL and N 7 LL) contributions of P +(4) ns in this limit. The corresponding coefficients read The C F n 4 f terms in eqs. (3.11) -(3.15) agree with refs. [94,95]; all other contributions have not been presented before. After transformation to x-space using eq. (2.23), the above results provide all small-x enhanced contributions to the non-singlet splitting functions P ± ns, L (x, α s ) in the large-n c limit at four and five loops. These results are shown in fig. 1 for an expansion in powers of α s , not a s = α s /(4π) as used in all formulae, at n f = 4 flavours. As seen before for the 3-loop splitting functions, see figs. 2 and 3 of ref. [5], and coefficient functions, see figs. 7 and 9 of ref. [8] -see also ref. [97] -all logarithms are needed for a meaningful approximation at any physically relevant small values of x. The N 7 LL small-x limit of P (4) ns, L (x) has already been employed in the first estimate of the 5-loop quark cusp anomalous dimension in the large-n c limit in ref. [96].
In this context, it is instructive to consider the generalization of the NNLL small-x resummation (3.5), which can be extended to N 7 LL accuracy in the large-n c (L) limit: x I 0 (z) + n 2 c (647 − 1008 ζ 2 + 144 ζ 3 + 1620 ζ 4 ) n c a 3 s ln 2 1 x I 1 (z) where we have omitted (with '. . .') the lengthy and (at least in the present notation) 'irregular' parts of the N 6 LL and N 7 LL coefficients for brevity; the ancillary file of this article provides the complete expression. The first two terms at each logarithmic order, which dominate at asymptotically small x for a fixed α s , can be seen to build up an exponential function, viz This exponential prefactor dampens, and from some unphysically large value of α s overwhelms, the small-x rise of the modified Bessel functions (3.7).
So far we have presented the 1/N n+1 ⇔ x 0 ln n x results arising from eqs. (2.16) and (2.21). We now turn to the double logarithmic 1/(N+p) n+1 ⇔ x p ln n x contributions which, as discussed in section 2, can be analyzed in a completely analogous manner at even p for P + ns and odd p for P − ns , and hence at all integer p for their common large-n c limit P ± ns,L .
Thus the complete Taylor expansions of the coefficients of ln k x can be determined at k = 6, 5, 4 for P (3) ns,L and at k = 2n, . . . , 2n − 3 for P (n) ns,L at n ≥ 4. In practice, we have performed the necessary computations to p = 70, which was more than sufficient to overconstrain the analytic expressions in terms of harmonic polylogarithms [83] for which an ansatz with up to 54 coefficients was used. This partial reconstruction of the analytic form of P ± ns,L can be carried out to any order in α s ; here we confine ourselves to the N 3 LO and N 4 LO expressions. The former is given by  (1 − x) −1 − 29 31 x H 1,1 p qq H 1,1 +

Results for the non-singlet coefficient functions
The order-by-order resummation of the highest 1/N powers of the even-N based ('+') non-singlet coefficient functions C + a for the structure functions in eq. (2.17) is 'automatically' included in the calculations towards the corresponding expressions for the splitting function P + ns reported in the previous section. Using the short-hand F ≡ S −1/2 with S as defined in eq. (3.1), the NNLL results for C + 2 and C + L analogous to eq. (3.2) for P + ns are found to be The corresponding result for the ν−ν charged-current (CC) structure function F 3 reads Like the splitting function, these coefficient functions exhibit a structure similar to that of their time-like counterparts in ref. [41], apart from the all-important fact that the sign of ξ is different there. The corresponding x-space results can be expressed in terms of generalized hypergeometric functions listed in appendix B. These hypergeometric functions related to the function F are, however, not present in the non-singlet physical evolution kernels K + a given by and . (4.10) The LL, NLL and NNLL x-space approximations resulting from eqs. (4.5), (4.6), (4.8) and (4.9) with eq. (2.23) are shown for four flavours in figs. 2 and 3. Unsurprisingly, these results alone do not provide relevant information about the behaviour of these coefficient functions at any physically interesting value of x. The results in eqs. (4.5) -(4.10) can become useful, however, once combined with other partial results such as the large-x limit and a sufficient number of moments.
As for the non-singlet splitting functions, the resummation can be performed for any even or odd power of x, and for all powers in the common large-n c limit C a,L of the coefficient functions C + a and C − a . Since the N 3 LO splitting functions are known in this limit [47], the all-x coefficients of the four highest small-x logarithms can be predicted at any higher order for a = 2 and a = 3.   (1 − x) −1 − 11 5 x + 48 5 x 2 − 72 5 x 3 H 0,1 − x . (4.12) We have also determined the corresponding five-loop results. Since they would become useful mainly in the context of research towards N 5 LO all-x expressions, which we do not expect in the foreseeable future, we skip these here for brevity but provide them in the ancillary file of this paper.

Results for flavour-singlet quantities
We now turn to the singlet case, and first present the results for the splitting functions P ik . The diagonal (i = k) quantities can be written as sums of 'non-singlet' (ns) and pure-singlet (ps) pieces, where P + gg is an non-singlet-like quantity, i.e., P gg in the limit C F = 0 (cf. ref. [98]). At leadinglogarithmic (LL) accuracy, Consequently the LL P + gg (x) is an oscillatory function, in notable contrast to the corresponding quark quantity in eqs. (3.2) and (3.5). The other LL contributions are found to be P ps (n) qq where the relation between P qg and P gq holds only at LL. In these expressions, ⌊. . .⌋ denotes the Gauß bracket (floor function), ρ = n − k − 2 i − 1, δ = n − k − 2 i, and C n are the Catalan numbers, (5.8) The first (second) binomial factors in eqs. (5.4) -(5.6) may be interpreted as counting the number of ways k (ρ/δ) gluons can be emitted from the gluon (quark) propagators, the number of which is directly related to the number of factors n f C F . In this interpretation, the quark and gluon emissions contribute to these logarithms equally, with emission strengths proportional to the colour factors.
Incidentally, the Catalan numbers (5.8) are the coefficients of the Taylor expansion of S in the LL non-singlet splitting function, suggesting the singlet results could be generalizations of their simpler non-singlet counterparts. However, we have not found closed-form expressions. The NLL and NNLL results are considerably more complicated, and even expressions in the form of the sum representations above are not available. To achieve this could require defining the LL results in terms of new special functions, and suitably generalizing them for the sub-leading results. Potential complications also arise from the observation that the NLL and NNLL results have terms with denominators similar to that of the N-space logarithmic functions in the SIA case [41], which could indicate the presence of logarithmic functions that depend on the special functions yet to be found.
The NNLL double-logarithmic resummation of the 1/N pole terms leads to the following results for the four-loop (N 3 LO) splitting functions (P +(3) ns (N) has been given in eq. (3.4) above) (1711 + 108 ζ 2 )C 3 A + 8 9 (2951 + 300 ζ 2 ) n f C F C A + 8 27 which we expect to become relevant in the near future in combination with the fixed-N moments in ref. [11] and other constraints. Their N 4 LO counterparts read P (4) qq (N) = P +(4) n 2 f C F + 32 3 (442 + 105 ζ 2 ) n f C 2 A + 32 9 (7054 (9211 − 72108 ζ 2 )C 2 F C A − 16 9 (16829 + 1602 ζ 2 ) n f C F C A − 8 9 (24337 − 22320 ζ 2 )C F C 2 [907 + 396 ζ 2 ]C 5 A − 640 9 [164 + 81 [3274 + 495 ζ 2 ] n 2 f C F C 2 A − 16 9 [38371 + 3978 Finally we present the corresponding results for the singlet structure functions F 2 and F L . The quark coefficient functions can be written as sum of non-singlet and pure-singlet pieces, The leading-logarithmic resummations of the pure-singlet and gluon coefficient functions are where ρ ′ = n − k − 2i − 2, δ ′ = n − k − 2i − 1, and the D n are defined as These are the coefficients of the Taylor expansion of i.e., those of the leading-logarithmic contributions to the non-singlet coefficient functions (4.1) - (4.3). Similarly to the splitting functions above, the first (second) binomial factor could be interpreted as the number of ways k (ρ ′ /δ ′ ) gluons can be emitted from the gluon (quark) propagators. Again, analytic NLL and NNLL results, which could require generalization of the LL singlet coefficient functions, are not available at this point. Corresponding results for the scalar-exchange structure function F φ in eq. (2.17), which is of only theoretical relevance, can be found in appendix C.
Hence, at least for the time being, the coefficients for the NLL and NNLL contributions to these coefficient functions are only known at each order separately. Since we do not expect yet higher orders to be become relevant in the foreseeable future, we finally present these results at the fourth and fifth order in a s = α s /(4 π). The former read (4913 + 11988 ζ 2 )C F C A + 16 27 (8461 − 1188 ζ 2 ) n f C F −   , and c (5) where the non-singlet contributions can be found in eqs. (4.8) and (4.9).

Summary and outlook
We have presented a comprehensive study of high-energy double logarithms appearing at the n-th order in perturbation theory in the splitting functions for the evolution of the parton distribution functions and in the coefficient functions of inclusive deep-inelastic scattering. These have the structure α n s x p ln 2n−n 0 −k x , where p ≥ 0 in gauge-boson-exchange DIS, and the parameter n 0 depends on the quantity under consideration. k denotes the logarithmic accuracy: k = 0 provides the leading-logarithmic (LL) terms, k = 1 the next-to-leading logarithmic (NLL) contributions etc.
For the flavour non-singlet quantities, the dominant contributions start at p = 0 with an offset of n 0 = 2 for the splitting functions P ± and the coefficient function C ± L , and of n 0 = 1 for the coefficient functions C ± 2 and C ± 3 . The structure of the unfactorized n-th order partonic structure functions in dimensional regularization has been employed to perform a NNLL resummation of these small-x logarithms to all orders in full QCD for the splitting function P + and the coefficient functions C + 2 , C + 3 and C + L ; the former can be expressed in terms of modified Bessel functions. Using, in addition, the known structure of the singular terms for P + in Mellin N-space as N → 0, the all-order resummation of the leading ln x terms has been pushed up to N 7 LL accuracy in the large-n c limit, where the functions P + and P − coincide. For the coefficient functions C + 2 and C + 3 the large-n c resummation has been extended to N 3 LL accuracy. In all cases, explicit fixed-order expansions up to the fifth order in perturbation theory have been presented for future reference.
In the flavour-singlet sector of standard DIS, the dominant small-x contributions at the n-th order are proportional to one inverse power of x enhanced by single logarithms, i.e., the splitting functions P ik and the singlet coefficient functions C 2 and C L at n loops are of the form x −1 ln m x with m = 1, . . . , n − n ′ 0 in the small-x limit. Our study has not added any new information on these terms or their all-order resummation, which is a long-standing and prominent, yet in its full generality still open problem in QCD. Instead, we have considered the double logarithms appearing with powers x p for even p ≥ 0, which correspond to expansions in Mellin N-space around (N + p) = 0. In these cases, we can again use dimensional regularization and exploit the structure of the unfactorized partonic structure functions at the n-th order in perturbation theory. We have presented NNLL small-x (p = 0) predictions up to five loops in full QCD. We can compute order-by-order to a very high power of α s in this framework, but we have not been able to find an all-order form which generates these results beyond the leading logarithms.
The results, certainly in the singlet sector, are not of immediate phenomenological relevance in DIS. However, they elucidate the analytic structure underlying the expressions at fixed order in perturbation theory. In addition, they provide important information for complete analytic computations of those quantities. In this regard, the results for the non-singlet splitting functions have already been used in the determination of the all-N expression of the large-n c non-singlet splitting function at four loops, based on a limited number of Mellin moments together with constraints on its endpoint behaviour and its the functional form. This application also allowed for important independent checks on the methods employed in the present article for the study of the high-energy (single and) double logarithms, and in turn provided input coefficients for the N 3 LL resummations.
While the intriguing structures of the resummations performed contribute to a much improved theoretical understanding, the chosen approach has also clear limitations. Most notably the leading p = 0 contributions to the non-singlet splitting function P − and the main ν+ν charged-current structure function F 3 are not accessible beyond the large-n c limit. Progress in this direction will require new methods in addition to those developed and considered here. The closed-form allorder resummation of the small-x double logarithms in the singlet sector remains an open problem, pending the identification of the proper set of functions, which complement the modified Bessel functions found to suffice in the non-singlet sector. Finally, the systematic study of DIS with an exchanged scalar and the implications for the flavour-singlet coefficient functions is a subject we have touched upon only briefly, with a few results presented in the appendix. These issues deserve further thorough investigation, which we leave for the future.
An ancillary files with our results in FORM format is available from http://arXiv.org.
Appendix A: expansions about N = 0 to order α 3 s Here we present the expansions of the fixed-order results used for the resummation of the x 0 ln k x contributions to the even-N based splitting and coefficient functions. The N n LL predictions are fixed by the corresponding N n LO results, hence only the N −4 coefficients of the fourth-order splitting functions P (3) are missing for the N 3 LL resummation of F 2 , F 3 and F φ . In view of a future determination of these quantities, the results below are given at N 3 LL accuracy.
The corresponding expressions for the LO, NLO and N 2 LO non-singlet splitting functions read The input coefficient functions for F 2 in eq. (2.7) in Laurent expansions analogous to eq. (2.22) are The corresponding expansion coefficients for the longitudinal structure function are given by and the input coefficient functions for the even-N based F 3 read Next we present the input quantities for the singlet cases. The 'diagonal' quantities can be separated into non-singlet and pure-singlet pieces. The quark cases in the following expressions are written in terms of the non-singlet parts presented above. Here we present the input expressions only to N 2 LL accuracy. The singlet splitting functions, expanded about N = 0 are given by C F C A − 12C 2 F + 16 3 C A n f + C A n f − + 16ζ 2 C 2 A + 76 9 C A n f − 32C F n f N −1 + O(N 0 ) , + 96ζ 2 C A C F n f − [120 − 32ζ 2 ]C 2 F n f − 16 9 C A n 2 f + The input coefficients c φ,i for F φ , which provides the resummation of P qg and P gg , are given by Appendix C: leading-logarithmic coefficient functions for F φ In this final appendix, we present some analytic results for the scalar-exchange coefficient functions C φ,q and C φ g . Using the decomposition C φ,g (N) = C + φ (N) +C where δ ′ and ρ ′ have been given below eq. (5.21).
Unlike the splitting functions and coefficient functions for F 2 and F L , the coefficient functions for φ-exchange DIS exhibit double logarithms also in BFKL-limit, i.e., they include contributions of the form α n s x −1 ln 2n−n 0 −k x . We have considered these in what might be called 'extended quantum gluodynamics' (eQGD), i.e., QCD in the limit C F = 0, and found an analogue to eq. (2.16) for the expansion about N = 1 which generates the resummation of these double logarithms.
In terms of S ′′ = (1 − 4 ξ ′′ ) −1 with ξ ′′ = C A a s /N 2 where N ≡ N − 1 , (C.5) the resummed coefficient function, including the finite NNLL contributions at O(a s ), is found to be where the two terms in the first line provide the LL and NLL contributions, and the remaining four lines the NNLL result. This resummation does, of course, also return the corresponding x −1 double-logarithms in P gg | C F =0 . These are found to vanish, as they have to, see refs. [51][52][53].