Small-x resummation from HELL

Small-x logarithmic enhancements arising from high-energy gluon emissions affect both the evolution of collinearly-factorized parton densities and partonic coefficient functions. With the higher collider energy reached by the LHC, the prospect of a future high-energy collider, and the recent deep-inelastic scattering (DIS) results at small-x from HERA, providing phenomenological tools for performing small-x resummation has become of great relevance. In this paper we discuss a framework to perform small-x resummation for both parton evolution and partonic coefficient functions and we describe its implementation in a computer code named High-Energy Large Logarithms (HELL). We present resummed and matched results for the DGLAP splitting functions and, as a proof of principle, for the massless structure functions in DIS. Furthermore, we discuss the uncertainty from subleading terms on our results.


Introduction
One aspect that makes the physics program of the CERN Large Hadron Collider (LHC) particularly rich is the vast kinematic region that can be explored. For inclusive enough processes, the kinematics is traditionally parametrized with a dimensionful scale Q, the typical hard scale of a process, e.g. a final-state invariant mass, and with the dimensionless ratio x = Q 2 /s, with √ s the machine energy. Thus, the success of the LHC physics program relies upon having control of the many ingredients that enter theoretical predictions, over a wide kinematic range in both x and Q 2 . This includes highorder corrections in QCD and in the electro-weak sector, resummation effects and non-perturbative inputs to hadronhadron cross section such as parton distribution functions (PDFs), which often represent the main source of theoretical uncertainty.
The bulk of experimental data that constrain PDFs comes from deep-inelastic scattering (DIS) data collected by the HERA experiments [1], which span several orders of magnitude in both x and Q 2 . Here, we concentrate on the highenergy, or small-x, regime. In particular, at low Q 2 , these data reach very small values of x, perhaps outside the region of validity of the fixed-order calculations which are used as inputs in the fits. Moreover, in the context of LHC physics, the unique design of the LHCb detector (essentially a forward spectrometer) makes this experiment well-suited to access a region of phase-space of very large rapidities, thus providing useful data to pin down the largely unconstraint PDFs at small x. The success of this enterprise relies on having a reliable theory description of the low-x region.
Despite the wealth of calculations listed above, very few phenomenological studies that incorporate both fixed-order and resummed calculations exist. The reason for this is technical: small-x resummation requires an all-order class of subleading corrections in order to lead to stable results. The purpose of this paper is to remedy this deficiency. We develop a framework to perform small-x resummed phenomenology. Our starting point is the resummation of coefficient and splitting functions according to the formalism developed by Altarelli, Ball and Forte (ABF) [17][18][19][20][21][22]. However, as we will describe in the paper, we introduce a number of improvements that make the procedure easier to extend to new processes, as well as numerically more stable. For the first time, we make resummed splitting and coefficient functions available in a public code named HELL (High-Energy Large Logarithms).
The structure of this paper is the following. In Sect. 2 we describe the ABF resummation of the splitting functions and its HELL implementation, highlighting and motivating several improvements. We then perform a comparison to the ABF original results and also to the ones of Ref. [16]. In Sect. 3 we introduce a method to perform the resummation of coefficient functions directly in transverse momentum space, which is then implemented in HELL. We show its equivalence to the ABF Mellin-space resummation, while discussing the numerous advantages of the new method. As a proof of principle, we present results for the partonic coefficient functions of the massless DIS structure functions F 2 and F L , as well as their comparison to the results obtained by ABF in Ref. [22].
Finally, we draw our conclusions in Sect. 4 and we outline forthcoming phenomenological studies which include fits of PDFs, as well as studies of small-x effects in electro-weak boson production at the LHC and Future Circular Colliders (FCC). Technical details are collected in a number of appendices.

Resummation of DGLAP evolution kernels
In this section we review the construction of resummed DGLAP evolution kernels needed for resummed PDF evolution up to NLL. We follow the formalism developed in the ABF series of papers [17][18][19][20][21][22]. We will also comment about other approaches, but leave a thorough analytic comparison to future work. Most of the section is devoted to introducing notation and describing how the theoretical results can be practically implemented in the code HELL. We will also present several improvements over the original implementation.
It is convenient to work in the space of the variable N conjugate by Mellin transformation to the variable x, since all convolutions become ordinary products. Here f i (x, Q 2 ) is a generic PDF, and we used a non-standard notation for the Mellin transform in which the kernel is x N rather than x N −1 . This is useful when discussing small-x because the small-x singularities, of the form (1/x) ln k x, are mapped into poles in N = 0 (in the usual notation, the poles are in N = 1): LL contributions at small-x correspond to terms in Eq. (2) with k = n to all orders n in α s , while NLL ones have k = n −1. Note that double logarithmic corrections, which would correspond to k = 2n, are absent in QCD, with the noticeable exception of the Higgs production in gluon fusion with a pointlike effective vertex in the large-m t effective theory [41].
The dominant small-x logarithmic enhancement only affects the singlet sector, while (double) logarithmic terms in the non-singlet are power-suppressed, i.e. they correspond to poles in N = −1. Therefore, we focus on the 2 × 2 singlet evolution matrix. The construction of resummed anomalous dimensions, which are the Mellin transform of the splitting functions, can be divided into three successive steps: 1. resummation of the "largest" eigenvalue γ + of the singlet anomalous dimension matrix 2. resummation of the quark-sector anomalous dimension γ qg 3. construction of the resummed anomalous dimension matrix in the physical (flavor) basis.
We address these three steps in turn, giving a brief summary of the ABF procedure, emphasizing those aspects that are different from the original construction. We finally comment on the numerical implementation and present some results.

Resummation of the largest eigenvalue
The singlet-sector DGLAP evolution equation reads where f g = f g (N , Q 2 ) and f q = f q (N , Q 2 ) are the gluon and quark-singlet PDFs respectively, and the evolution matrix is given by (omitting arguments for readability) As already mentioned, the non-singlet sector is not affected by small-x logarithmic enhancement, and we therefore ignore it. The DGLAP evolution equation Eq. (3) can be diagonalised by performing a change of basis. We define the "eigenvectors" f ± as where the transformation matrix R (and its inverse) can be generically written as Substituting Eq. (5) into Eq. (3) we get In general, to make the equation diagonal, one has to provide a matrix R such that the matrix in squared brackets in Eq. (7) is diagonal, Solving this problem in general is rather complicated. However, we notice that at pure LL level the matrix that diagonalizes has constant coefficients, so we can ignore the second term in squared brackets and simply solve an eigenvalue problem. At NLL, a non-trivial dependence on Q 2 appears; however, the action of the derivative with respect to Q 2 further suppresses the second term in squared brackets by α s β 0 , showing that it first contributes at NNLL level. Therefore, when treating running coupling effects perturbatively, we can ignore the derivative contribution and simply focus on the eigenvalue problem, which in particular leads to the following explicit form for R, being γ ± the eigenvalues of . We anticipate that running coupling effects will eventually be resummed to all orders in α s β 0 : when this counting is adopted, the derivative term is no longer subleading and the matrix R should be corrected for it. We will come back to this point later in Sects. 2.3 and 3.2.
The eigenvalue γ + is chosen to be the largest eigenvalue at small-x, i.e. N ∼ 0, namely the one which is enhanced at small N , while γ − is finite in N = 0. Consequently, f + is the only eigenvector that contains logarithmic enhancement and which is affected by high-energy resummation. This holds for several factorization schemes, including DIS and MS, and the so-called Q 0 MS scheme which is particularly useful in small-x resummation [32,34,46,47]. The resummation of small-x logarithms in the evolution is then encoded in the resummation of the largest eigenvalue γ + . The difference between the MS and Q 0 MS factorization schemes influences the resummation of γ + beyond the leading logarithmic accuracy, as well as the resummation of γ qg and of the coefficient functions, as we shall see in more detail in Sect. 3. The structure of the resummation described in the remainder of the section is rather general and it is valid for both MS and Q 0 MS schemes. When presenting phenomenological results our scheme of choice will be Q 0 MS, which is preferred from an all-order viewpoint, because it gives more stable results [22]. It has to be noted that, when expanded to fixed-order, the difference between the two schemes only starts at relative O(α 3 s ): thus, all theoretical predictions that enter current PDF fits are not sensitive to this choice.
High-energy resummation is achieved thanks to the BFKL equation [7][8][9][10][11][12], which, in analogy with DGLAP, we write as an evolution equation for the moments of the parton density. Therefore, defining the M moments of f + by with Q 0 some reference scale (the PDFs depend logarithmically on Q, so the value of Q 0 is irrelevant), we have where χ is the BFKL kernel, currently known to NLO [12] and to NNLO in the collinear approximation [47] (see Ref. [48] for recent work beyond NLO accuracy). In the small-x and high-Q 2 limit, both the DGLAP and BFKL equations are expected to hold, and consistency between the solutions to both equations allows to resum to all orders collinear contributions in the BFKL kernerl or, equivalently, small-x contributions in the DGLAP anomalous dimension.
Knowledge of the BFKL kernel to N k LO accuracy allows for the resummation of the N k LL contributions to the DGLAP anomalous dimension (and vice-versa). It is worth noting that Eq. (11) is an ordinary differential equation only if the coupling does not run. Indeed, in M-space, α s (Q 2 ) becomes a differential operatorα s , essentially because ln Q 2 is turned into −∂/∂ M and consequently Eq. (11) is to be intended as an operator-valued equation. This is a manifestation of the well-known fact that the eigenvalues of the LO kernel do not diagonalize the BFKL equation at NLO. Consistency between DGLAP and BFKL equations allows us to build a double-leading (DL) expansion of γ + and χ which takes into account the logarithmically enhanced contributions in both ln Q 2 and ln(1/x) [13]. Because of the poor perturbative behaviour of the BFKL kernel, obtaining a stable resummed result is however not straightforward and requires a somehow complex procedure with a careful treatment of the formally subleading terms. This issue received great attention in the past, mainly by three groups: Refs. [13][14][15][16], Refs. [17][18][19][20][21][22][23], and Refs. [24][25][26][27]. Despite the different approaches, which are characterized by different treatments of formally subleading corrections, a fairly consistent picture emerged, with small differences in the final results of the different groups (see e.g. Ref. [49]).
The ABF approach [17][18][19][20][21][22], which we adopt in this paper with a few improvements, allows us to build perturbatively stable resummed results by combining four main ingredients: duality, symmetrization, momentum conservation and running coupling resummation, as we summarize below.
Duality between the DGLAP anomalous dimensions and the BFKL evolution kernel, is the statement that in the fixed coupling limit (i.e. neglecting contributions due to the running of α s ), the kernels satisfy the following relations [50,51] Beyond LL Eq. (12) is corrected by contributions due to the running of α s . In principle Eq. (12) provides all the ingredients for small-x resummation: we start with the BFKL kernel χ at a given order (LO or NLO) and we use duality to determine a DGLAP anomalous dimension, dual to χ , which resums small-x contributions to the desired logarithmic accuracy (LL or NLL). However, as previously mentioned, the BFKL kernel itself exhibits a very poor perturbative behaviour, with poles of the form α k s /( j − M) k for any integer j at every perturbative order k. The poles in M = 0 and M = 1, which correspond to the collinear and anti-collinear regions, are particularly harmful [13]. The key observation is that the resummation of collinear poles (which in momentum space are just collinear logarithms) is controlled by the DGLAP anomalous dimension. Hence, we can use duality, in the opposite direction, to derive a kernel χ , dual to standard DGLAP, that resums all the collinear enhancements. The DL kernel can then be constructed by matching standard BFKL with the collinearly improved one. Furthermore, again by duality, this result can be turned into an anomalous dimension.
However, the stabilization of the collinear region does not completely cure the problem, because of the singularity of the BFKL kernel in M = 1. Indeed the behavior in middle region between M = 0 and M = 1 determines by duality the nature of the rightmost small-N singularity, i.e. the asymptotic small-x behaviour of the splitting functions. The nature of the singularity obtained in this way is perturbatively unstable: it is a pole at fixed order, a square root branch-cut at DL-LO, non-singular at DL-NLO, see e.g. [52]. The anticollinear terms can however be resummed and thus stabilized by exploiting the symmetry properties of the BFKL kernel, which relate them to the collinear contributions [13,21]. This symmetrization is performed by constructing a kernel which coincides with the DL one at a given logarithmic accuracy in ln Q 2 and ln(1/x), but satisfies the required symmetry properties exactly (while in general these would be spoiled by subleading terms). In the ABF approach, the symmetrized kernel is defined via implicit equations which must be solved numerically (more details are given in Appendix A). Note that the definition of the symmetrized kernels has some degree of arbitrariness, due to the inclusion of unconstrained subleading contributions. After symmetrization, the singular behaviour of the dual DGLAP anomalous dimension is always a square root branch-cut.
The third important ingredient of the ABF resummation is momentum conservation, which implies that the first Mellin moment of the largest eigenvalue must vanish, and translates by duality into a constraint on the BFKL kernel: In general, in a DL expansion, Eq. (13) is violated by subleading terms, but it may be enforced by adding a subleading contribution which does not introduce new singularities at small N and vanishes at large N . We stress that, while the final anomalous dimension must clearly satisfy momentum conservation, one can decide whether momentum conservation should be imposed in any of the intermediate steps.
We note that the stability of the result is greatly improved by enforcing momentum conservation at each step of the resummation procedure. Symmetrization and momentum conservation allow us to build perturbatively stable BFKL kernels, and, by duality, DGLAP anomalous dimensions, in the fixed coupling limit. The resulting singularity is however modified at every perturbative order by running coupling corrections to duality [19]. These corrections start at NLL and, while formally subleading, they are in fact dominant since they change the nature of the small-N singularity. The dominant running coupling corrections are resummed by solving the running-coupling BFKL evolution equation for f + , and extracting its anomalous dimension (see e.g. Refs. [24][25][26]). This can be done analytically by approximating the kernel in proximity of its minimum, which in turn corresponds by duality to the squareroot branch cut of the anomalous dimension, i.e. its leading singularity. After running coupling resummation, the rightmost singularity of the anomalous dimension is turned back to a simple pole (as it was at fixed leading order), but now shifted from N = 0 to N = N B (α s ) > 0. The overall effect is a suppression of the small-x growth with respect to the (symmetrized) DL result.
Combing all the effects together, the final form of the resummed DGLAP eigenvalue in the ABF approach at LO+LL is while at NLO+NLL is In the above equations γ ,(N)LO + contains the symmetrized double-leading contributions at LO and NLO respectively, which include the fixed-order part of the anomalous dimensions. The "Bateman" contribution γ B,(N)LL contains the running coupling effects obtained by solving the evolution equation, and carries the actual small-N singularity. The remaining term in each equation avoids double counting. Further details and explicit formulas are given in Appendix A. For later convenience, we also define which contain only the resummed contributions to be added to the corresponding fixed order. Note that one could also imagine to match the resummation to NNLO. This step, which is usually straightforward, is rather cumbersome in this case essentially because the dependence on the strong coupling of symmetrized DL result γ + is only known numerically. We leave this further matching for future work, stressing that it is of great interest especially in the context of PDF fits. Before moving to the resummation of the quark and gluon entries of the anomalous dimension matrix, let us briefly comment about the different approaches to the resummation of γ + that can be found in the literature. The resummation proposed in Refs. [13][14][15][16] is based on very similar ingredients as ABF, namely the resummation of collinear singularities, symmetrization and running coupling effect. However, rather than relying upon duality to determine the resummed anomalous dimension, the running coupling BFKL equation is solved and the anomalous dimension is extracted from the solution. References [24][25][26][27] on the other hand, only relied on running coupling corrections and not on symmetrization, with the argument that high-Q 2 physics should be dominated by the M ∼ 0 region. A thorough study of all the sources of uncertainty in small-x resummation of γ + would require investigating all the subleading modifications described above and goes beyond the scope of this work. However, given the conclusions of Ref. [49], which found the three different approaches to be in reasonable agreement, one might expect the small-x resummation of the eigenvalue γ + to be under good control.

Resummation of the quark anomalous dimension
The high-energy behaviour of the qg anomalous dimension has been derived at the leading logarithmic level in Refs. [33,34]. The quark anomalous dimensions are always suppressed by a power of α s with respect to the gluon ones, so they enter for the first time at NLL.
The all-order small-N behaviour of γ qg is determined from the resummed anomalous dimension γ + : In order to perform the resummation, the function h to all orders in its argument is needed; however, to the best of our knowledge, a closed form for h in either MS or Q 0 MS does not exist. Nevertheless, the coefficients h k of its Taylor expansion can be computed recursively, as described in Ref. [34]. The first 35 coefficients have been worked out in Ref. [22]. The singular behaviour of γ qg up to O(α k s ) is obtained by including the singular behaviour of γ + up to at least the same order.
We first address the question of which accuracy is needed for γ + in Eq. (17). Since NLL effects in γ + will contribute to NNLL in γ qg , we could use the LL expression for the largest eigenvalue. However, since the position of the pole determines the asymptotic small-x behaviour of the result, the use of the LL γ + pole is not ideal because it would lead to displaced poles in different entries of the anomalous dimension matrix. Therefore, we find it convenient (mostly from a numerical point of view) to use an hybrid expression which we denote LL which is based on the DL-LO result but contains the running-coupling NLL contribution. In formulae, we define In other words, this expression is basically the same as γ LO+LL + , Eq. (14), but the parameters entering the Bateman anomalous dimension γ B (and consequently all the double counting terms), which determine the position of the pole, are those of the NLL result Eq. (15).
The function γ LO+LL + , Eq. (19), cannot be directly used in Eq. (17), because its growth at large N (due to its fixed-order component) would produce a spurious large N behavior in γ qg to all orders in α s . Therefore, we use where γ LO,sing + is the singular N ∼ 0 part of the LO anomalous dimension. 1 We point out that this procedure differs from that of Ref. [22], where h is computed with γ NLO+NLL + , and the large-N behaviour is subtracted by recomputing h with γ NLO + − γ NLO,sing + . We comment on the differences between the two approaches in Appendix B.2. Here we just stress that the two procedures are formally equivalent, our formulation leading to a faster and more reliable numerical implementation.
The resummation of running coupling contributions also affects the determination of γ qg . In the approach of Ref. [35], it is included by computing 1 In principle it would be sufficient to include in only the singular LL contributions. However, one might argue that it is safer to also include additional subleading (NLL) terms in it, provided they vanish at large N . We indeed include these NLL terms; details are given in Appendix B.2.
where the square brackets notation γ k is defined by the recursion the dot denoting the derivative with respect to ln Q 2 . In our implementation,γ is computed as a derivative with respect to α s ,γ = −β 0 α 2 s ∂γ /∂α s . The need to compute a derivative with respect to α s of the resummed anomalous dimension is one of the main practical motivation for using γ LL + rather than γ NLL + , as the numerical evaluation of the former is much faster and more stable than the latter, thereby allowing a more precise determination of the numerical derivative. Note that Eq. (22) comes from the approximate assumption that γ is linear in α s [35] (we will explicitly re-derive this result in the context of coefficient functions in in Sect. 3.3). Under this assumption,γ would simply beγ −β 0 α s γ . We shall also consider this additional approximate form forγ as a means to estimate the uncertainty due to this approximation.
A further complication arises from the fact that after the inclusion of running coupling corrections Eq. (21), the series Eq. (18) is only asymptotic. In Ref. [22] the resummation is performed by computing the sum of the series à la Borel, using a truncated Borel integral corrected with an asymptotic behaviour derived from a simpler solvable model. We adopt here a different approximate procedure, which only relies on the available information from h. We make use of a Borel-Padé summation procedure, where we compute the sum of the series à la Borel, and use a Padé approximant for the sum of the Borel-transformed series obtained from a finite number of coefficients of the expansion of h. Details of this procedure are given in Appendix B.1.
Finally, from Eq. (21) we can construct the pure resummed contribution as the contribution to be added to the NLO anomalous dimension to obtain a matched NLO+NLL result.

Construction of the resummed singlet splitting function matrix
Now that we have resummed the largest eigenvalue and the qg component, we can construct the full anomalous dimension matrix in the gluon-singlet basis. First of all, the qq component can be recovered by making use of the colorcharge relation [34] γ NLL The eigenvalue γ − , which is finite in N = 0 and does not resum, contains a finite fixed-order constant terms which is formally NLL, This particular form, together with the color-charge relation Eq. (24), is such that the r − component of the transformation matrix, Eq. (9), is simply given at NLL by The gg component can be recovered by transforming back the diagonal matrix to the physical basis, leading to the general expression Using Eq. (9), valid in the fixed-coupling case, we simply get γ gg = γ + + γ − − γ qq , which combined with Eqs. (24) and (25) leads to When resumming running coupling effects, the form of r ± changes and, consequently, Eq. (27) receives in principle running-coupling corrections. However, we have checked that these effects are typically smaller than the various sources of ambiguity in the whole resummation procedure coming from subleading contributions. Therefore, following Ref. [22], and without loss of accuracy, we adopt Eq. (27) as our default implementation for γ gg . On the other hand, a more careful treatment of running coupling effects is needed when dealing with the resummation of coefficient functions, as we shall see later in Sect. 3.2. Finally, it remains to compute γ gq ; however, the available information is not sufficient to constraint its NLL part. This is not a problem, because the accuracy of the solution of the evolution equation is formally NLL even if the gq entry is just LL. At LL, we can just use a color-charge relation this equation can be modified by using the NLL expression of γ gg , even though the resulting gq anomalous dimension will still remain formally accurate at LL. For phenomenological application we find useful to write the resummed and matched anomalous dimensions as a fixedorder contribution γ (N)LO plus a γ (N)LL , which contains the resummation minus double counting. Thus, in this notation, the NLO+NLL evolution matrix is given by where γ NLL qg is given in Eq. (23), and γ NLL gg = γ NLL can be easily derived from Eq. (27). From the above matrix, one can compute the inverse Mellin transform and obtain the resummed splitting functions, where P NLL gg and P NLL qg are the ultimate primary ingredients for a resummed DGLAP evolution.
The results in momentum space deserve further comments. The contributions γ i j vanish, by construction, at large N . This is enough to guarantee that their x-space conjugates P i j are ordinary functions, i.e. they do not contain plus distribution or delta functions. However, they potentially exhibit a constant behavior, or even an integrable singularity, as x → 1. To avoid potential problems with matching at large x, we follow ABF and we further suppress these functions at x = 1 with an x-space damping: However, despite the many desirable features of the above damping procedure, momentum is no longer conserved in Eq. (29). In the flavor basis, momentum conservation implies that namely, the sum of each column must vanish in N = 1. Both equations imply which is violated. The origin of the violation is twofold: first, even though γ NLL + is originally constructed to vanish in N = 1, it looses this property once the x-space damping is applied; second, γ NLL qg is not necessarily vanishing at N = 1 even in absence of damping. While this momentum violation was not considered in the original ABF work [22], here we force momentum conservation by a modification of the gg entry and d(N ) is any function that goes to zero at large N and has no leading singularities. We use so that the momentum conservation can be restored directly in both N and x space. Given the numerous steps involved in the resummation procedure, we find useful to summarize our strategy in implementing them in a numerical code: -we compute γ NLL The four P i j constructed in this way are the primary output of the code HELL.

Numerical implementation and results
The numerical implementation of the resummation of γ + is quite challenging. The main difficulty comes from the fact that several ingredients of the resummation procedure are not available in a closed analytic form, but they are only defined as zeroes of implicit equations which must be solved numerically in the complex plane. Moreover, these equations can depend on functions which are themselves computed as zeros of implicit equations (see Appendix A for more details and explicit examples). While for real N one can rely on robust root-finding algorithms such as bracketing methods, in the complex plane one must rely on root-polishing methods whose convergence heavily depends on the accuracy of the initial guess supplied to the algorithm. Moreover, several functions have more than a single branch which satisfy the zero criterium, hence it is crucial to consistently identify the correct one.
We circumvent the above difficulties by computing γ (N)LO+(N)LL + (N , α s ) only along the contour for Mellin inversion, which we parametrize, in the upper plane Im N > 0 (in the lower plane we use the complex conjugate path), as is the integration variable and c ∼ 1 is a parameter whose value is adjusted for each value of α s to give optimal convergence properties for the Mellin inversion. For t = 0, N = c is real, and we can therefore use robust bracketing root-finiding algorithms which are guaranteed to converge. As we move from N = c into the complex plane (t > 0), we resort to the secant method, whose reliability entirely depends on our ability to provide an accurate guess of the root to be found. Our strategy here consists in proceeding by small steps in t, using for initial guess at each step the value of the function at the previous step. If the step is fine enough and the function sufficiently well behaved, this method works well and also avoids jumps across different branches. Very rarely, when this method fails, we can also use a slower but more stable minimum-finding algorithm, by turning the problem of finding a zero of a function into the one of finding the minimum of the absolute value of the function itself. As a consistency check, we verify that at large |N | (large t) the resummed expression becomes asymptotically close to the known fixed-order result.
Using this strategy, we construct tables of values of γ (N)LL + (N , α s ) along the contour for a grid in α s , one grid for each value of n f = 3, 4, 5, 6. The tables also contain information about the leading singularities of γ + , namely the position of the leading poles and value of their residues. We keep the code which produces the tables private, and use the tables as primary ingredients for the public code presented in this work.
The public code HELL reads the provided tables as input files, and performs the remaining steps for the resummation. In particular, it constructs the resummed quark anomalous dimension γ NLL qg (N , α s ) according to the procedure described in Sect. 2.2, along the Mellin inversion contour. It then performs the inverse Mellin transform and reconstruct the full singlet splitting function matrix, as described in Sect. 2.3. (A similar strategy is used for the resummed coefficient functions, see Sect. 3.) The HELL code, while being quite flexible and numerically stable, is rather "heavy" (∼ 100 MB) due to the size of the files which contain the tabulated γ (N)LL + , and also slow due to the presence of numerical integration (although we implemented a dynamical caching which speeds up multiple evaluation in a single run). Therefore, we created a higherlevel variant of the code, dubbed HELL-x, which reads pretabulated (with HELL) splitting functions (and coefficient HELL-x has been interfaced to the evolution code APFEL [53], and will be in future used to obtain high-energy resummed PDF fits. We now present some representative results for the resummed splitting functions for α s = 0.2 and n f = 4. In Fig. 1 we show all four entries of the evolution matrix: P gg (upper-left panel) and P gq (upper-right panel), P qg (lowerleft panel) and P qq (lower-right panel). The values of x range from 1 to 10 −9 . We include in the plots the fixed-order splitting functions at LO (dashed), NLO (solid) and NNLO (dot-dot-dashed) in black. At resummed level, we show in solid purple the NLO+NLL result, while the LO+LL result is shown in dashed green and is present only for P gg and P gq , as the other two entries do not have any leading logarithmic enhancement. At NLO+NLL we have to specify the factorization scheme. As previously mentioned, we adopt Q 0 MS, which is convenient from an all-order viewpoint. We recall that the difference between MS and Q 0 MS starts relative order O(α 3 s ) and therefore the fixed-order splitting functions start to differ only beyond NNLO.
We see that at large x the resummation has no effect, due to the damping, so the resummed result smoothly matches onto the fixed order. At smaller x, the resummed result grows. The effect is more pronounced in the case of P qg , where the growth starts immediately, while for P gg the growth is delayed by an initial decrease, a well-known feature of subleading small-x contributions [15,21,27].
Similarly, we see the same effect on P qq and P gq , where the contribution of the resummation is just C F /C A times the contribution on the left plots, Eq. (30). As far as P gg and P gq are concerned, we observe a nice perturbative convergence of the resummed and matched results, with the NLO+NLL being a very small correction over the LO+LL, especially when compared with the fixed-order perturbative behaviour at small x. This convergence derives from the stability of the resummation of γ + , mostly determined by the the constraints imposed by symmetrization and momentum conservation, as described in Sect. 2.1.
We have included in Fig. 1 an "uncertainty band" for the NLO+NLL result. This band is determined by replacing in Eq. (22)γ with −α s β 0 γ . As Eq. (22) is derived under the assumption of linearity of γ LL + , both expressions are equally valid, and the difference between the two can be taken as a measure of the uncertainty coming from subleading corrections beyond the linear approximation. The distance between our default construction and this alternative approach is then symmetrized, thus giving the band. We acknowledge that the resulting uncertainty is just one of the many sources of uncertainties of the resummation, as coming from the various approximations described before and from subleading terms. However, we think that the uncertainty shown in Fig. 1 is a good representative of the uncertainty from subleading contributions. We have indeed verified that other variations of subleading terms, e.g., the actual form of γ LO,sing + in Eq. (20), leads to similar effects. On the other hand, the uncertainty on the resummation of γ + is likely to be much smaller, due to the many constraints on its construction, as confirmed the agreement between different groups [49], as well as by the good convergence of the gluon entries. Clearly, the overall uncertainty from all sources of ambiguities will be larger, but we believe the shape and the relative size among the various entries is likely to be well represented by the current band.
We now move to the comparison of our results with other approaches. To better highlight the impact of the resummation, we show the comparisons in terms of ratios over the fixed-order splitting functions. In Fig. 2 the ratio of resummed LO+LL splitting functions over the LO ones are presented for P gg and P gq (at this order, only the gluon components are affected by resummation). Along with our curves, the ABF results of Ref. [22] are also shown in dotdashed cyan (the plotted range is limited in x due to the available information from the original paper). The fixed NLO (solid) and NNLO (dot-dot-dashed) are also shown (in gray) for comparison's sake. Overall, we observe good agreement with our result. The tiny deviation is due to a different treatment of the n f dependence of the result, see Appendix A for more detail. Interestingly, we observe that at large x the resummed results tend to follow the shape of the NLO and NNLO results, before merging onto the LO due to the damping, perhaps an indication that higher order contributions predicted by the resummation go in the right direction even far from the small-x region. Note also that the LO+LL ratio is basically identical for P gg and P gq , a small difference being visible only at large x. This is easily understood by noting that the small-x behaviour of both fixed-order and resummed results are simply related by a color factor C F /C A .
The comparison of the NLO+NLL resummed results are shown in Fig. 3. Here, not only we compare our results to the ones obtained by ABF in Ref. [22] but also to the resummed splitting function calculated in Ref. [16] (henceforth the CCSS approach). The latter also comes with a (yellow) uncertainty band which is obtained from renormalization scale variation. While the agreement with ABF is still rather good, there are more significant deviations, especially in the quark entries, which come from many sources. For P qg (and P qq ), we use the LL anomalous dimension, Eq. (20), while ABF used the full NLL anomalous dimension. Moreover, we implement differently the large-N subtraction, as discussed in Sect. 2.2, and we also have different numerical implementations, as we adopt a Borel-Padé summation for the series Eq. (21). These differences also affect P gg (and P gq ), due to Eq. (27) but their numerical impact appears to be smaller. Note that for these gluon splitting functions we also have differences at large x due to our implementation of . Unfortunately, our simple uncertainty band does not fully cover all these differences, especially at larger x. When comparing to CCSS, we see that the gluon entries P gg and P gq are in decent agreement, our result lying at the lower edge of the CCSS band. The quark entries P qg and P qq , however, are quite different both in shape and in size. It is clear that these entries are affected by larger uncertainties, as demonstrated by both our and the CCSS bands, as well as by the large perturbative corrections in the fixed order. Therefore, it is likely that such a difference is a manifestation of this ambiguity, which could be fixed only by a NNLL computation. Note also that the CCSS results are obtained in a scheme which is not exactly the Q 0 MS, and it is well known that differences between schemes can be significant at resummed level (see e.g. the comparison of the MS and Q 0 MS in Ref. [22]).

Resummation of perturbative coefficient functions
We now turn our attention to the resummation of small-x enhanced contribution to collinearly factorized partonic coefficient functions. The general formalism for the resummation of inclusive cross sections is based on k t -factorization, which was derived a long time ago [29][30][31][32][33][34] and it is known to LL 2 for an increasing number of cross sections and distributions [36][37][38]40,41,[43][44][45]54].
The ABF approach for resumming coefficient functions was developed in Ref. [35] and applied to the case of DIS structure functions in Ref. [22]. The crucial point to note is that, analogously to the case of PDF evolution, the resummation of formally subleading running coupling corrections plays a crucial role. The procedure that we will describe in this section does take these effects into account but departs from the original ABF method in that the resummation is performed directly in transverse momentum space rather than in Mellin moment space. Although the two procedures are formally equivalent, as we shall discuss below, the momentumspace technique significantly helps with two shortcomings of the Mellin-space approach. First of all, computing Mellin moments of k t -factorized cross sections with respect of the gluons' k t often constitutes the bottle-neck of a calculation. Secondly, running coupling corrections in Mellin space are included order-by-order in perturbation theory and then a Borel summation of the resulting series is performed, resulting in potential numerical instabilities. Thus, working directly in transverse-momentum space avoids dealing with asymptotic series and opens up the possibility of performing resummed calculations for processes for which Mellin moments cannot be computed analytically.
In order to keep the notation simple, we consider a process with only one hadron in the initial state, such as DIS. The generalization to two hadronic legs is straightforward, as discussed in Ref. [38]. Because we are interested in the high-energy limit, we limit ourselves to consider the singlet sector. Although we work in transverse-momentum space, we find convenient to take Mellin moments with respect the longitudinal momentum fractions and to work with cross sections in N space. The generic cross section is then given by where f q and f g are the quark-singlet and gluon PDFs, respectively, and in the second line we have transformed to the basis of the eigenvectors of singlet DGLAP evolution. In DIS, σ can be either the structure function F 2 or F L (F 3 is non-singlet), up to a normalization factor (for precise definitions, see Ref. [34]). Since only f + (N , Q 2 ) resums at small x, we have a single coefficient function, C + (N , α s ), which is affected by small x enhancements. We will come back later in Sect. 3.2 on the precise definition of C + in terms of C g and C q . It is known, e.g. [29][30][31], that in the high-energy limit a different, more general, form of factorization holds, even away from the collinear limit 3 : is the unintegrated (k t dependent) gluon PDF and C(N , k 2 t /Q 2 , α s ) is the off-shell coefficient function, i.e. the coefficient function for the partonic process with an off-shell initial state gluon.
In the high-energy limit, the unintegrated gluon density can be related to the standard resummed PDF Before discussing the form of U(N , k 2 t /Q 2 ), we immediately observe that once the relation Eq. (39) between the integrated and unintegrated PDFs is established, by comparing the gluon contribution in k t -factorization Eq. (38) and the high-energy contribution in collinear factorization Eq. (37) we are able to write This equation represents our main formula for the implementation of high-energy resummation in the coefficient functions. In HELL, the k 2 t integral is evaluated numerically, given the off-shell cross section C(N , k 2 t /Q 2 , α s ) in k t space, and an actual form of U(N , k 2 t /Q 2 ), which will be discussed in the next subsection. Note that LL accuracy only requires to calculate C to lowest order in α s . Moreover, its N dependence is also subleading and one can set N = 0.

The evolution factor
We now turn to discussing the form of U(N , k 2 t /Q 2 ) in Eq. (39). As clear from Eq. (39), it first evolves the largest eigenvector PDF from Q 2 to the scale k 2 t , where it then converts it to the unintegrated gluon PDFs. It can be understood either in terms of the all-order gluon Green's function [29][30][31]34] or as the evolution kernel of a generalized ladder expansion [36].
At lowest order and fixed coupling, the form of U is known [34] U s N , where γ s is the anomalous dimension obtained from the leading order BFKL kernel with duality at fixed coupling α s = α s (Q 2 ). We also note the scheme-dependent factor R(γ s ) that originates from the correct treatment of collinear singularities, the calculation of which requires a more accurate analysis away from d = 4 space-time dimensions. In the commonly used MS scheme this factor reads [34] where χ 0 (M) is the eigenvalue of the leading-order BFKL kernel and (x) and ψ(x) are the Euler gamma and digamma functions, respectively. In Q 0 MS instead we simply have Comparing the last line of Eq. (42) to Eq. (43) we see that the difference between the two schemes starts at relative O α 3 s . It is also useful to write the scheme-dependent factor as while it is not straightforward to find a closed analytic form ofR MS (ξ ) from Eq. (42), in the Q 0 MS scheme we simply The running-coupling generalization of Eq (41) that we implement is where γ + is the resummed anomalous dimension. Note that by substituting γ + → γ s , at fixed coupling, we recover the lowest-order result Eq. (41). The general structure of the result appears fairly complicated because of the presence of the scheme factorR. However, in the preferred scheme Q 0 MS the first two lines of Eq. (45) evaluate to unity and the result simplifies to where we recognize the derivative of a DGLAP evolution factor.

Basis transformation and collinear subtraction
Once C + is computed according to Eq. (40), one can use the relation between C + and C g , C q to obtain resummed expressions for C g and C q . This relation can be trivially obtained from the transformation matrix that diagonalizes the DGLAP evolution equation in the singlet sector, Eq. (6), leading to At fixed coupling, as discussed in Sect. 2.1, diagonalizing the evolution equation simply amounts to diagonalizing the singlet anomalous dimension matrix, and using Eq. (9) would lead to the simple relations However, as previously discussed, when running coupling effects are taken into account, a transformation that diagonalizes the evolution matrix does not in general diagonalize the evolution equation, since the derivative with respect to Q 2 acts on the transformation matrix, Eq. (7), producing an additional contribution which is in general not diagonal. Furthermore, we note that in contrast to the case of γ gg , here a more careful treatment of these running coupling corrections is required in order to guarantee the all-order cancellation of collinear singularities that may be present in C + . Finding the general transformation matrix that diagonalizes the singlet evolution equation is not an easy task. However, because our goal is to find a running coupling version of Eq. (48), a full solution is not needed, as long as we limit ourselves to the LL accuracy.
To this purpose, it is convenient to consider the logarithmic derivative of Eq. (37) with respect to Q 2 The first two and last two lines of Eq. (49) are related by the same transformation matrix that relates first and second line of Eq. (37). However, the logarithmic derivative already produces running coupling contributions, making further running-coupling effects on the transformation matrix genuinely subleading. Thanks to this observation, we can use the fixed-coupling transformation matrix to relate the various terms in Eq. (49). For the + component we are mostly interested into, this leads to the equation We now need to understand the logarithmic order of each contribution, and keep only those terms which are LL. First, we observe that the logarithmic derivative of the coefficient function is one logarithmic order higher than the coefficient function itself. This suggest that all derivative terms could be thrown away, leading back Eq. (37). However, the key point of the resummation of running coupling effects is exactly to keep those subleading terms which are suppressed by α s β 0 , which are precisely those coming from these derivatives. Next, from the analysis of the previous section, we know that γ gg and γ gq are LL, while γ qg and γ qq are NLL.
Since to all orders C q is of the same logarithmic order as C g (as we shall see later in this section), this suggests that only the first two terms on the right-hand side of Eq. (49) should be kept. However, some of those terms can be leading if there is a fixed-order contribution in the coefficient function which is of higher logarithmic order than the coefficient function itself. This is for instance the case of the DIS structure function F 2 : in this case, both C g and C q are NLL (in absolute counting), but the fixed-order expansion of C q is where the first O(α 0 s ) term is formally LL. When this is the case, the term C q γ qg with C q replaced by its fixed-order superleading contribution leads to a leading contribution to the equation and must be retained. Finally, the last contribution is genuinely subleading.
After all these consideration, and further approximating γ gg with γ + (the difference being subleading), we end up with the equation which can be easily solved introducing an exponential factor so that Eq. (51) becomes The solution is then having defined where Q 0 is the scale at which U vanishes (which is the position of the Landau pole), and we have left C q outside the integral because it is either 1 or 0. Eq. (54) represents the running coupling version of the first of Eq. (48), at LL. As a cross check, we can easily verify that if the coupling does not run we get which is indeed equivalent to Eq. (48) up to subleading terms. With similar arguments, it is also possible to show that the solution in presence of running of the equation for C − leads exactly to its fixed-coupling counterpart, second line of Eq. (48), up to NLL terms. Note that this suggests that the generalization of the transformation matrix R, Eq. (6), is simply obtained by using (up to subleading corrections) and the fixed-coupling value of r − .
We have now all the ingredients to obtain resummed expressions for C g and C q . From Eq. (54) we immediately have where in the second line we have used Eq. (40). As we already discussed, the C q subtraction is suppressed by a NLL term, so this term is present only when C q has a fixed-order contribution which is superleading. This is the case of the DIS structure function F 2 , where C 2,q is NLL (in absolute counting) and C 2,q = 1 + O(α s ). In this case, we can just replace C q with 1 and get In other cases, such as the longitudinal structure functions F L , C L ,q is still NLL in absolute counting but does not contain any superleading fixed-order contributions, as it starts at O(α s ); therefore, the C q contribution is genuinely subleading and one finds The resummed expressions for C q can be found from the second of Eq. (48), We note that C − does not contain any logarithmic enhancements to all orders and therefore it can be safely evaluated at fixed-order (NLO) and at N = 0. Furthermore, we can make use of the high-energy color-charge relation Eq. (24) and arrive at which shows that C q and C g are of the same logarithmic order, as anticipated.
In order to perform the matching to the fixed-order, we find useful to introduce (i = g, q) where C (k) i (N ) is the kth order coefficient of the α s expansion of the resummed result C i (N , α s ). Hence, the second contribution subtracts from the first term (the resummed result) its expansion up to the perturbative order we want to match to, e.g. n = 1 is NLO, n = 2 is NNLO. In this notation the resummed and matched contribution is simply given by where the resummed contributions n C i (N , α s ) are computed by HELL, while the fixed-order parts have to be provided by an external code. Note that the color-charge relation Eq. (62) reduces to when written in terms of n contributions, provided n ≥ 1. Note also that these n C q , with n ≥ 1, can be seen as the resummed contributions to the pure-singlet quark coefficient functions [34,38].

Equivalence between transverse-momentum space and Mellin space resummations
In this section we want to compare our result Eq. (40) with running coupling Eq. (45) to the analogous result obtained in the ABF approach [22], which is performed in Mellin space. For convenience, and without loss of generality, we work in the Q 0 MS scheme, and using the definition Eq. (52) we can write where we introduced the dimensionless variable ξ = k 2 t /Q 2 . We can thus write Eq. (40) as In the ABF approach the resummation of coefficient functions closely follows the one of the quark anomalous dimension γ qg , where in place of the function h(M), Eq. (18), the Mellin transform of the off-shell coefficient function with respect to k t is used. Therefore we define the so-called impact factor, 4 whereC(N , M, α s ) admits an expansion in powers of M Note that k ≥ −1 for processes that are not two-particle irreducible in the high-energy limit and therefore their lowestorder diagrams exhibit a collinear singularity, as in the case of F 2 , while k ≥ 0 for processes without such collinear singularity, as in the case of F L . The inverse Mellin transform is given by where the integration contour is a standard Mellin inversion contour, with 0 < c < 1. The resummed expression for the coefficient function C + can be now found substituting Eq. (70) into Eq. (67). The integral over ξ can be performed in all cases and we find where we have introduced a lower integration limit ξ 0 . This lower limit is equal to 0 in the fixed coupling case, but in the running coupling case we have due to the presence of the Landau pole. Note that, assuming γ + > 0 (as appropriate close to the pole), U (N , ξ 0 ) = 0 so Eq. (71) simplifies which represents an equivalent form of Eq. (40).
Let us now focus on the simpler case without collinear singularities,C −1 = 0. We want to show that the sum in Eq. (73) corresponds to the procedure adopted in ABF, under some assumptions on the form of the resummed anomalous dimension. In particular, we recover ABF assuming that the dominant running coupling effects are determined by 1-loop running of the lowest power of α s appearing in the anomalous dimension. In other words, one makes the approximation [as which is an exact expression at LO, where γ + N , α s (μ 2 ) = α s (μ 2 )γ (0) + (N ). In order to better describe the exact anomalous dimension which is not simply linear in α s , one can replace so that the μ 2 derivative of Eq. (74) in μ 2 = Q 2 is correct (and the 1-loop structure is kept). In this particular approximation, the ν-derivatives in Eq. (73) satisfy the recursion where U ABF indicates the evolution factor Eq. (52) computed with γ + from Eq. (74). We recognize the recursion defined in Eq. (22). This recursive construction is exactly the method employed by ABF to perform the running coupling resummation of coefficient functions. Therefore, we recover the ABF result 5 (in the case of no collinear singularities, as in F L ) where we further computed the expansion coefficientsC k in N = 0. However, we recall that the resulting series is divergent, and cannot be summed analytically, so sophisticated numerical techniques with limited numerical accuracy are needed in order to use Eq. (73), see Appendix B.1.
In presence of a collinear singularity, the first term in Eq. (73) proportional toC −1 (N , α s ) does not vanish. Additionally, the collinear subtraction due to C q must be included. In the ABF approach, the subtraction is written first in Mellin space as α s h(M)/M, with h(M) defined in Sect. 2.2, and subtracted directly at the level of inverse Mellin integrand, leading to (in the case of F 2 ) To prove the equivalence of Eq.  79) it is immediate to verify that the h 0 term cancels, and the integral can be computed to lead to exactly Eq. (78). Note that the usage of the running-coupling version of the basis transformation discussed in Sect. 3.2 is crucial to obtain the correct result. Had one used the fixedcoupling version, the collinear singularity would not cancel. Therefore, we have shown that our transverse momentum space derivation and the Mellin space resummation adopted by ABF are completely equivalent, even though the current result is more general and does not rely on the assumptions of Eqs. (74) and (75). A numerical comparison of Mellin-space and k t -space resummation is performed in Fig. 4. The plot shows 2 C a,g Eq. (63) for both a = 2, L, with α s = 0.2 and n f = 4 in the Q 0 MS scheme. We observe that the two approaches give indeed the same result. We note however that the Mellin space implementation suffers from numerical instabilities, which determine small oscillations around the actual result. These oscillations become more severe at larger α s , and disappear at smaller α s . We note that these numerical instabilities are related to the approximate Borel-Padé method used for the Mellin space implementation, which necessarily uses a limited amount of information (i.e., a finite number of coefficients of the M = 0 expansion, see Appendix B.1). In Ref. [22] a different "truncated" Borel method was used, which did not develop oscillation; however, also in that case the amount of information used was limited, while in our k t -space approach we make use of all the information residing in the off-shell cross section.

Numerical implementation and results
We now turn to the numerical implementation of the resummation of coefficient functions in HELL. Starting from Eq. (40) written as in Eq. (67), we integrate by parts (the boundary terms vanish at ξ → ∞ thanks to C and in ξ 0 thanks to U ) and evaluate the off-shell cross section at N = 0 (since its N dependence is subleading), As the resummation of coefficient functions is at present accurate only at LL, we may conveniently compute U (N , ξ) using the LL anomalous dimension introduced in Eq. (20), However, since α s in the evolution factor is evaluated at Q 2 ζ with ζ ranging up to ξ , and ξ is integrated over all accessible values, the resummed anomalous dimension should be computed at extreme values of α s , from 0 to ∞. This is problematic in practice, since the resummed anomalous dimension is itself computed numerically as described in Sect. 2, and it is numerically challenging to reach both high and low values of α s . Therefore, a convenient implementation consists in adopting the approximation Eq. ( This expression is advantageous because the integral in the evolution factor has been computed analytically and it only requires γ LL + and its α s derivative at α s = α s (Q 2 ). We now turn to the specific case of massless DIS. For an observable without collinear singularity, such as the longitudinal structure function, we simply have For processes with collinear singularities, we further need the collinear subtraction U qg , Eq. (55), to obtain C g , Eq. (58).
Computing the integral Eq. (55) numerically, even within the approximation Eq. (74), is challenging due to the need of integrating γ qg over a range of α s from α s (Q 2 ) to ∞. In principle, we could find an approximation similar to Eq. (74) for The sum in Eq. (86) can be computed as we compute γ qg itself. In fact, the computation is identical, except that the h k coefficients are all shifted by a unity. This way, we can pre-tabulate it once for all, and use it for any observable with collinear singularities. The integral term in Eq. (86) can be combined with the integral in Eq. (81), so that the collinear subtraction is performed at the level of the integrand, leading to a more reliable numerical implementation. So, for C 2 , we have finally From these resummed expressions, we can then construct the resummed contributions, n C g (N , α s ), Eq. (63) (see also Appendix C.2), and n C q (N , α s ) from Eq. (65). At this point, as we did for the splitting functions, we damp the resummed contributions in x space multiplying by (1 − x) 2 to ensure a smooth matching onto the fixed order. 6 6 In practice, a smoother matching to NLO is obtained if 1 is derived from 2 , as detailed in Appendix C.
The resummed and matched partonic coefficient functions are shown in Fig. 5 in the case of C L , and in Fig. 6 in the case of C 2 . In both cases, the gluonic coefficient functions are shown on the left-hand panel, while the quark ones on the right-hand panel. The solid purple line is for NLO+LL, while the solid green for NNLO+LL. The resummation is performed in Q 0 MS. Analogously to the case of the splitting functions, the size of the uncertainty band is obtained from the symmetrized difference between the calculation performed with r as given in Eq. (84) or its linearized version r = α s β 0 . The corresponding fixed-order results are also shown: NLO in dashed, NNLO in dot-dotdashed and N 3 LO [56] in dotted. The plots are for α s = 0.2 and n f = 4.
The comparison to the ABF approach is done in Fig. 7, where the resummed contribution 1 C a,g (a = 2, L) is shown. We note that our results are in general agreement with the ones of the ABF paper [22], especially if we focus on the longitudinal coefficient functions C L ,i , i = g, q. In the case of C 2,i , differences are instead more pronounced. This should not come as a surprise because, as discussed at length, the resummation for the coefficient functions differs by various subleading terms. We stress once again that we have verified (see e.g. Fig. 4) that the resummation performed in Mellin space (as in Ref. [22]) gives identical results (modulo numerical instabilities at large α s ) as our k t -space formulation, as long as the same γ + is used and the same subtraction of the large-N terms is adopted. Therefore, the difference comes from both the different way of subtracting the large-N behaviour (see discussion in Appendix B.2) and the fact that we use γ LL + rather than γ NLL + . Moreover, note that the band is indeed larger in the C 2 case, confirming that subleading effects in C 2 are more pronounced than in C L . In particular, Fig. 6 Same as Fig. 5, but for the coefficient functions C 2,i Fig. 7 Comparison of the resummation contribution 1 C a,g as obtained in this work (solid) versus the ABF results (dot-dashed) of Ref. [22] for both a = 2 (purple) and a = L (green), using α s = 0.2 and n f = 4 in the Q 0 MS scheme a direct comparison with the expressions of Ref. [22] shows that our result differs by constant terms at O(α s ) and O(α 2 s ) in the resummed γ + , which lead to formally NLL and NNLL differences in the resummed coefficient functions. We conclude that, in absence of a well motivated preference for these subleading contributions, both results have to be considered as equally valid at the present logarithmic accuracy, and the ambiguity can only be fixed by computing (resumming) the NLL contributions in the coefficient functions. At larger x, we observe a significant deviation between our result and ABF for C L , which is not well represented by the band. In this case the difference has to do with the large-N matching, and we expect our matching procedure to perform better than ABF.

Conclusions and outlook
In this paper we have discussed the resummation of highenergy, i.e. small-x, logarithms that affect both the evolution of collinearly-factorized parton densities and perturbative coefficient functions. Despite a wealth of calculations have been performed in k t -factorization, the framework that allows for high-energy resummation, very few phenomenological studies that incorporate both fixed-order and resummed calculations existed, essentially because of the complexity of the running-coupling resummation of the DGLAP and BFKL evolution kernels.
In this paper we have overcome this obstacle and we have developed a computer code named HELL (High Energy Large Logarithms), available for download at www.ge.infn.it/~bonvini/hell , that enables one to obtain small-x resummed DGLAP splitting and partonic coefficient functions. The code is based on the formalism developed by Altarelli, Ball and Forte (ABF), with several improvements that avoid numerical instabilities and facilitate the future inclusion of different processes. The main innovation with respect to the ABF original procedure consists in performing the resummation of perturbative coefficient functions from the off-shell cross section in transverse-momentum space rather than in Mellin-moment conjugate space. Therefore, partonic off-shell cross sections computed in k t -factorization can be directly used, without the necessity of performing Mellin transformations with respect to the initial-state gluons' virtualities, which is often the bottle-neck of this kind of calculations.
We have provided resummed results for the splitting functions in the singlet sector, both at LO+LL and NLO+NLL and, as a proof of principle, we have also performed the resummation for the massless DIS structure functions F 2 and F L , at NLO+LL and NNLO+LL. We have provided a qualitative estimate of the theoretical uncertainty by varying subleading contributions that are related to the running of the strong coupling. We have found that this uncertainty is rather small for the gluon splitting functions P gg and P gq , essentially because their resummation is mostly driven by the all-order behaviour of the leading eigenvalue in the singlet sector, which is under good theoretical control. On the other hand, the uncertainty is larger for the quark splitting functions P qg and P qq , as well as for the closely-related DIS coefficient functions, for which we only control the first tower of logarithmic contributions. This feature also appears in the comparisons to ABF and CCSS. Indeed, all the approaches considered here are in decent agreement for the gluon splitting functions, while they significantly differ in the quark sector, which is also plagued by rather large uncertainties.
We see this, rather technical, paper as the first necessary step towards a rich program of small-x phenomenology. First, we would like to use the results presented here to perform a PDF fit of DIS data that consistently include small-x resummation in both parton evolution and perturbative coefficient functions, especially in view of the recent final release of HERA data [1]. These small-x resummed PDF fits will be performed in the NNPDF global analysis framework [57] and preliminary results have been presented in [58].
Furthermore, having at hand resummed PDFs, we will perform a study of small-x effects at high-energy hadron colliders, such as the LHC or an FCC. In particular at FCC, because of the extremely large center-of-mass energy, low-x effects in processes like Higgs or vector boson production are expected to become very important. In this respect, the study of electro-weak boson production via the Drell-Yan mechanism offers an almost unique environment to look for deviation from standard DGLAP dynamics. Finally, some of us have recently developed frameworks to combine small-x resummation with threshold [59] and transverse-momentum resummation [60] and we look forward to performing phenomenological studies of joint resummation in the context of

A Details on the resummation of γ +
In this appendix we give further details about the resummation of γ + . A comprehensive treatment of this topic can be found in Refs. [21,22] and it is beyond the purpose of this paper. Here, we collect the relevant formulae needed for the numerical implementation of our version of the ABF resummation procedure, pointing out the changes and improvements we introduced with respect to the literature.

A.1 Double-leading contributions and symmetrization
As we briefly explained in Sect. 2.1, one of the ingredients for building a stable DL expansion of the BFKL kernel (and by duality of the resummed DGLAP anomalous dimension) is symmetrization [13], i.e. the construction of a kernel which satisfies symmetry properties otherwise spoiled by subleading terms. As explained e.g. in Refs. [13,21], the kernel χ in the fixed coupling limit satisfies χ(M, α s ) = χ(1− M, α s ) if the kinematic is symmetric upon exchange of the virtualities Q 2 and k 2 . This is e.g. true for gluon-gluon scattering where the kinematic is x = Q 2 k 2 /s, but the symmetry is broken for DIS-like kinematics where x ∼ Q 2 /s. If χ σ and χ are the kernels obtained with a symmetric (gluon-gluon scattering) and asymmetric (DIS) choice of x respectively, one can however show that the following (equivalent) relations hold In the ABF approach one constructs, at a given logarithmic accuracy, a symmetric kernel χ σ such that and a corresponding asymmetric kernel χ satisfying Eq. (88), by means of the introduction of so-called off-shell kernels. An off-shell kernel is a kernel χ(M, N , α s ) which depends on both M and N and is related to its on-shell counterpart χ(M, α s ) and to the dual anomalous dimension γ (N , α s ) by the equation The actual expressions for the kernels are more involved and take into account several technical details which are thoroughly discussed in Ref. [21]. In the following we collect explicit formulae for the results. The anomalous dimension γ ,LO (14) is obtained by putting on-shell the kernel where here and in the following χ mom is a subleading contribution which implements momentum conservation. This can be taken to be of the form where f (0) = f (∞) = 0 and f m (1) = 1, e.g.
and c m is such that the final kernel satisfies momentum conservation In Eq. (96), χ s represents the dual of the LO DGLAP anomalous dimension and it is defined by the implicit equation The kernelχ 0 (M, N ) contains the BFKL LL contributions and is constructed from the LO BFKL kernel by performing the off-shell extension discussed above and removing the double counting with the DGLAP contributions. Its expression reads where the function γ s is the dual of the LO BFKL kernel The anomalous dimension γ ,NLO,fc + is obtained by putting on-shell the kernel where In the last equation χ s,NLO is the dual of the NLO DGLAP anomalous dimension Our construction of the kernelχ 1 (M, N ) follows the procedure outlined in Ref. [21] (see in particular Appendix A of that reference), but the result slightly differs since we separate the collinear and anti-collinear singularities in the whole range −∞ < M < +∞ rather than just on a finite interval (the impact of this on the resummed splitting functions is however very small and formally of higher twist). The result can be written as where and The symmetrized kernel χ 1 is written in terms of the function φ + L defined by as In the numerical implementation of the resummation procedure for n f = 0, in Eqs. (100) and (109) we do not use the exact eigenvalue of the DGLAP matrix γ + . The reason for this is the presence of a branch-cut singularity in the solution for the eigenvalue equations. Although this branch-cut cancels out in results for physical observables, in practice the exact cancellation is spoiled by subleading terms introduced in the resummation procedure. Since the cut is on the right of the leading small-N singularity, it introduces an unphysical oscillating behaviour in the splitting functions. One can however observe that the whole resummation procedure can be consistently carried out by replacing γ + with any function which has the same small-N behaviour. In our approach, we replace γ + with the same function evaluated in n f = 0 plus the n f dependence of the LL and NLL contributions only. As usual, we also add a subleading term which enforces exact momentum conservation. The only missing pieces from the resulting DL expansion are thus the n f -dependent contributions which are not enhanced at small N , but these always cancel out when taking the difference between resummed and unresummed result. The final result for γ LL + and γ NLL + defined by Eq. (16) is thus correct at the given logarithmic accuracy in both ln(1/x) and ln Q 2 but free of spurious branch cuts. In Ref. [22] a slightly different method was used, where a rational approximation of the whole n f dependent part was used, but we find our minimal approach cleaner and more convenient (note e.g. that by adding too many terms to the approximation one reconstructs an approximation of the branch cut and thus reproduces the unphysical oscillating behaviour of the result). We verified that the difference between the two approaches is numerically very small (and of course formally subleading). A minor subtlety arises when dealing with (subleading but de facto dominant) running coupling effects, which we discuss in the next section.

A.2 Running-coupling contributions
The leading small-N singularity of the anomalous dimension is determined by running coupling corrections which, as already mentioned, determine the small-x asymptotic behaviour of the splitting functions and can be resummed by solving the BFKL evolution equation for f + [24][25][26]. In the ABF approach [21], the resummation of the dominant running coupling contributions is encoded in the so-called Bateman anomalous dimension γ B,(N)LL (N , α s ) appearing in Eqs. (14) and (15) (henceforth generically referred to as γ B (N , α s )). The latter is determined from the solution of the BFKL evolution equation for f + , obtained from a quadratic approximation of the BFKL kernel around its minimum M = M min (α s ), which in turn corresponds to the rightmost singularity of the DL anomalous dimension. This implies that γ B (N , α s ) depends on the intercept c(α s ) and curvature κ(α s ) of the kernel in M = M min and their derivatives with respect to α s . These are referred to as Bateman parameters. The Bateman parameters are most conveniently computed using the BFKL kernel in symmetric variables, which is related to the kernel in DIS variables by Eq. (88).
Notice that the values of c and κ are the same for both kernels. The actual BFKL kernel used for the calculation of the Bateman parameters differs from the one discussed in the previous subsection in two respects. As we explained in The expressions in Eqs. (119) and (120) are new, since in Ref. [21] they were given as Taylor expansions around M = 1/2 (notice that M min − 1/2 = O(α s )). It is worth pointing out that these expressions, on top of removing the truncation error present in the mentioned Taylor expansion, can be easily evaluated numerically since multiple derivatives of χ s can always be written in terms of derivatives of γ LO + evaluated at the solution for χ s .
The result for the Bateman anomalous dimension can be written as where U (a, b, z) is the confluent hypergeometric function of the second kind and Notice that the only difference between γ B,LL and γ B,NLL is the kernel used for the computation of the Bateman parameters, as explained above, while their functional form is identical. An equivalent representation in terms of Bateman functions (from which the name for γ B ) is given in Ref. [21]. The double-counting terms between the Bateman anomalous dimension and the DL expansion which appear in Eqs. (14) and (15) can be written as at LO and NLO respectively, where The γ match term in Eqs. (127) and (128) removes a subleading spurious branch-cut due to using different kernels for the DL and the Bateman anomalous dimensions. It can be chosen to be of the form where c and κ are the Bateman parameters while c β 0 and κ β 0 are the intercept and curvature in the minimum of the final off-shell kernels defined in Eqs. (96) and (105) for the symmetrized DL result. Finally, the term γ mom is a subleading contribution which enforces momentum conservation and can be chosen to be of the same form as χ mom in Eq. (97). We finally observe that, for the anomalous dimension γ LO+LL + defined in Eq. (19), the double counting term γ LO,NLL d.c. has the same form of γ LO, d.c. but with the Bateman parameters c and κ computed from the NLO Bateman kernel (118) in order to match the singularities of γ B,NLL .

B Details on the resummation of γ q g
In this section we provide some detail on the resummation of γ qg . Note that what follows also applies to the resummation of U qg , Eq. (86), and of the partonic coefficient functions.

B.1 Borel-Padé method
Our starting point is either Eqs. (17) or (21), both of which provide the resummation of γ qg in terms of the function h(M) which is not known in closed form. Only the first few coefficients of its Taylor expansion in M, Eq. (18), are known. However, the usage of a truncated series will inevitably decrease the all-order logarithmic accuracy to a fixed-order accuracy. Therefore, we need a method to keep the all-order nature of the result, while dealing with just a finite set of coefficients.
The idea used here (originally proposed in Ref. [62]) is to construct a Padé approximant of the sum of the series from a given number of coefficients. In practice, given that the series is expected to diverge [22], the actual implementation consists in summing the series à la Borel, using a Padé approximant for the Borel-transformed series. Namely, Eq. (21) becomes where the inner sum is to be replaced by its Padé approximant. In Ref. [62] a diagonal [ p/ p] Padé, in which the degree of the numerator is identical to the degree of the denominator, was used. Here, we have found that a better numerical stability is achieved by using an almost-diagonal approximant [ p/( p − 1)], where the degree of the denominator is lower by a unity. We also consider a stronger second-order Borel summation, which basically consists in applying the Borel method twice, leading to [22,62] γ NLL where K 0 is a modified Bessel function of second kind. Again, the inner sum is to be replaced by its Padé approximant. While this stronger method might not be strictly necessary, we find that it performs better than the first-order method. Therefore, we use Eq. (135) for all applications presented in this work. This Borel-Padé method, though far from optimal, works reasonably well, provided α s is not too large (α s 0.3) and the number of coefficient used is not too high (we use p = 8, i.e. 16 coefficients). We adopt this method also for the computation of the part of U qg , Eq. (86), which is given as a series.

B.2 Large-N subtraction
We now briefly comment on the two ways of subtracting the large-N behaviour in the computation of γ qg we mentioned in Sect. 2.2. For ease of notation, let us denote such that γ LL + = γ −γ , and lim N →∞ Ignoring for simplicity the complications coming from the resummation of running coupling effects, which is not relevant for the present discussion, we implement the resummation of γ qg as Equation (21), which automatically vanishes in the large-N limit due to Eq. (137) (except for the 0-th order term of the series, which is anyway subtracted when matching to fixed order). This differs from the choice performed by ABF [22] where the first term contains the small-x resummation, and the second term subtracts the large-N behaviour by recomputing h withγ as argument; finally, the zero-th order constant is restored with the last term. 8 Note that in the original ABF formulation the full NLO+NLL anomalous dimension was used, however here we are interested in the different ways of subtracting the large-N behaviour, so we do not need to add this complication to the discussion. In other words, we apply the large-N subtraction before acting with h, while ABF do it after h. Our option is closer to the "plain" resummation obtained by γ s qg = h(γ s ), and has the advantage of having to compute h (through the Borel-Padé method) only once.
To understand the differences between the two approaches we expand the two results from which we can write the difference as γ ABF qg − γ qg =γ (γ −γ ) 2h 2 + 3h 3 γ + . . .
These terms vanish at large N because of Eq. (137), as they should, so the large-x limit is the same with the two procedures. At small N , close to the pole, these terms vanish only if γ vanishes in N = 0. When this is the case, it is then clear that the two approaches will give equivalent results (this is indeed what we find). However, ifγ does not vanish in N = 0, the difference, though formally subleading, can be sizeable. From this observation it seems favourable to constructγ such that it vanishes in N = 0. This is achieved if γ LO,sing + contains not only LL terms (as formally strictly necessary) but also NLL terms (which are not formally needed for the present accuracy). Expanding the LO at small N up to NLL we have γ LO,sing + = α s π Note that in the second term (the NLL contribution) we have added a damping factor 1/(N +1). This is needed because this NLL term is originally a constant, and therefore if included without damping it would spoil Eq. (137), and hence the large-N subtraction.

C Details of DIS resummation
In this appendix we collect the relevant expressions for the massless off-shell coefficient functions in DIS, needed for the resummation described in Sect. 3, and discuss the matching to fixed order.

C.1 Massless off-shell coefficient functions
The off-shell cross section in the case of massive quarks has been computed in Ref [30]; more precisely, the N = 0 moment of the cross section is Eq (4.12). 9 Here we take the massless limit of the expression reported in Ref. [30], obtaining C 2 (0, ξ, α s ) = n f α s 3π Its Mellin transform is [34] For the longitudinal coefficient function, we find an expression similar to Eq. (144) C L (0, ξ, α s ) = n f α s 3π 9 N moments are computed with respect to ρ = 4m 2 s rather than z = Q 2 s , however the difference is subleading.
its Mellin transform reproduces the result of Ref. [34]: For our numerical implementation, we find useful to note that both C 2 (0, ξ, α s ) and C L (0, ξ, α s ) satisfy where we made explicit use of the fact that C (N , ξ, α s ) is linear in α s , and we let ξ 0 → 0 since we are expanding to fixed order. The expansion of C 2 , Having the expansion of the resummed coefficient functions up to O(α 2 s ), we can now construct both 1 C i and 2 C i in N space, and then in x space by Mellin inversion. We have however noted that, while 2 C i (x, α s ), which contains contributions starting at O(α 3 s ), vanishes fast enough at large x (after applying the (1 − x) 2 damping discussed in Sect. 3.4) and hence ensures a smooth matching onto the fixed order, 1 C i (x, α s ), which contains contributions starting at O(α 2 s ), is sizable at large x ∼ 10 −1 , thereby potentially spoiling the accuracy of the resummed and matched NLO+LL result in that region. Since the culprit of this sizeable effect is exactly the O(α 2 s ) term of the expanded resummation, we find it convenient to re-define 1 C i (x, α s ) as where the damping (1 − x) 2 is the standard damping adopted everywhere, and f (x) is a further damping function such that f (0) = 1 and f (1) = 0. We have identified a convenient form for the damping function in f (x) = (1 − √ x) 6 . We observe that for f (x) = 1 one would recover the original definition.