Towards parton distribution functions with small-$x$ resummation: HELL 2.0

In global fits of parton distribution functions (PDFs) a large fraction of data points, mostly from the HERA collider, lies in a region of $x$ and $Q^2$ that is sensitive to small-$x$ logarithmic enhancements. Thus, the proper theoretical description of these data requires the inclusion of small-$x$ resummation. In this work we provide all the necessary ingredients to perform a PDF fit to deep-inelastic scattering (DIS) data which includes small-$x$ resummation in the evolution of PDFs and in the computation of DIS structure functions. To this purpose, not only we include the resummation of DIS massless structure functions, but we also consider the production of a massive final state (e.g. a charm quark), and the consistent resummation of mass collinear logarithms through the implementation of a variable flavour number scheme at small $x$. As a result, we perform the small-$x$ resummation of the matching conditions in PDF evolution at heavy flavour thresholds. The resummed results are accurate at next-to-leading logarithmic (NLL) accuracy and matched, for the first time, to next-to-next-to-leading order (NNLO). Furthermore, we improve on our previous work by considering a novel all-order treatment of running coupling contributions. These results, which are implemented in a new release of HELL, version 2.0, will allow to fit PDFs from DIS data at the highest possible theoretical accuracy, NNLO+NLL, thus providing an important step forward towards precision determination of PDFs and consequently precision phenomenology at the LHC and beyond.


Introduction
The outstanding quality of the CERN Large Hadron Collider (LHC) data continuously challenges the particle physics theoretical community to perform refined calculations with uncertainties comparable to the ones of the experimental results. Consequently, perturbative predictions for LHC processes nowadays often include radiative corrections in QCD at next-to-next-to-leading order (NNLO) and, in some cases, N 3 LO [1,2]. These remarkable calculations have a tremendous impact on many aspects of LHC phenomenology. This can be seen, for instance, in the context of the determination of parton distribution functions (PDFs), where the inclusion in the fit of data points describing the transverse momentum of Drell-Yan lepton pairs [3] or various differential distributions in top pair production [4] has recently become possible thanks to the existence of fully differential calculations at NNLO. Moreover, the recent completion of the NNLO QCD corrections to jet production [5,6], paves the way for global determinations of parton densities which are truly NNLO. On the other hand, it is well known that fixed-order calculations fail to provide reliable results in regions of phase space which are characterized by the presence of two or more disparate energy scales. In such cases, large logarithmic contributions appear to any order in perturbation theory, and they must be accounted for to all orders. Resummed calculations have also seen remarkable progress in recent years, and for many observables the resummation of next-to-next-to-leading (NNLL) logarithms of soft-collinear origin has become standard, even reaching N 3 LL for a few observables [7][8][9][10][11][12][13][14][15][16].
Traditionally, resummation is not included in the calculations that are employed in PDF fits with the argument that the observables which are considered are rather inclusive and therefore not much sensitive to logarithmic enhancements. However, the LHC is exploring a vast kinematic region in both the momentum transferred Q 2 and the Bjorken variable x = Q 2 /S, where √ S is the centre-of-mass energy. It is therefore important to assess the role of logarithmic corrections both in the small-x, i.e. high-energy, regime and at large x, i.e. in the threshold region. For instance, the production of a lepton pair via the Drell-Yan mechanism, which is measured by the LHCb collaboration in the forward region and at low values of the leptons' invariant mass, probes values of x down to 10 −5 ÷ 10 −6 . At the other end of the spectrum, searches for new resonances at high mass are sensitive to PDFs in the x ∼ 10 −1 region.
In the past, some of us included threshold resummation in PDF fits [17] (see also Ref. [18]) and performed dedicated studies that included threshold resummation in both coefficient functions and PDFs in the context of the production of heavy supersymmetric particles [19]. It has to be noted that the inclusion of threshold resummation in PDF fits did not pose particular challenges because in the widely used MS scheme the splitting functions that govern DGLAP evolution are not enhanced at large x [20,21], and the effect of threshold resummation can be included through a K-factor. The situation is radically different if we consider small-x resummation, where both coefficient functions and splitting functions receive single-logarithmic corrections to all orders in perturbation theory.
In this paper we further improve on our recent work, with the goal of providing all the necessary theoretical ingredients and numerical tools to perform a PDF fit which includes small-x resummation in both PDF evolution and in DIS partonic coefficient functions, including the correct treatment of the transitions at the heavy flavour thresholds, at NLL accuracy matched, for the first time, to NNLO. This is important because DIS data still represent the backbone of any PDF fits, and HERA data [62] do explore the small-x region. Achieving this target requires two ingredients which are the main results of this work. The first is the construction of a so-called "variable flavour number scheme" at the resummed level [63], which provides coefficient functions with powerbehaving mass effects and resummation of mass collinear logarithms, as well as describing the transition of PDF evolution at heavy quark thresholds. The second is the matching of the NLL resummed splitting functions to their fixed-order counterparts, which is realized for the first time up to NNLO. This requires the expansion in power of α s of the resummed result, which proved to be non-trivial because of the way the resummed kernels are constructed. We also correct an error present in the our previous paper [61] that we inherited from the original ABF work [37], which has however a small phenomenological impact. More importantly, we derive a new way to perform the resummation of running coupling effects, which streamlines the construction of the resummed evolution kernels and solves an issue in the construction of its α s expansion. All these improvements are included in a new release of HELL, version 2.0, which is publicly available at the webpage www.ge.infn.it/∼bonvini/hell.
We point out that the investigation of the impact of the resummation on the NNLO fixed-order splitting functions is potentially of great phenomenological interest. It is well known that, due to accidental zeros, the effect of small-x logarithms in the evolution is mild at NLO but stronger at NNLO, and will be even stronger (with two extra logarithms) at N 3 LO. On top of this, comparison of theoretical predictions with experimental data suggests that fixed-order NLO theory describes the small-x region better than NNLO, see e.g. [64]. Indeed, the final resummed anomalous dimensions tend to have a shape which is somewhat closer to NLO than to NNLO, as the small-x growth is greatly reduced once resummation is performed, see e.g. [65]. Therefore, having the possibility of using NNLO theory stabilized at small x with high-energy resummation can provide a significant improvement in the description of the data. This is indeed observed in preliminary applications of our results [66]. Furthermore, a reliable theory of DIS at low x finds interesting applications beyond collider physics. For instance, it is a key ingredient in the description of ultra high-energy neutrinos from cosmic rays, see e.g. [67][68][69].
This paper is organized as follows. In Sect. 2 we discuss the implementation of a variable flavour number scheme in the context of small-x resummation, while the details of resummation of DGLAP evolution and its expansion are collected in the two following sections. Specifically, in Sect. 3 we present our new realization of the resummation of running coupling contributions, and later in Sect. 4 we perform the perturbative expansion of the various ingredients and construct our final resummed and matched splitting functions. We present our numerical results in Sect. 5, before concluding in Sect. 6. In Appendix A we provide analytical results for the off-shell DIS partonic cross section with mass dependence, while in Appendix B we collect some technical details on the actual implementation of the resummation in HELL.

Resummation of deep-inelastic scattering structure functions
The standard way of describing the deep-inelastic scattering of an electron off a proton is to express the cross section in terms of structure functions that depend on the Bjorken variable x and the momentum transfer Q 2 , where sum runs over all active flavours, i.e. all the partons for which we consider a parton density. The number of active partons depends on the choice of factorization scheme and will be discussed in Sect. 2.1. Note that the coefficient functions also depend on the charges of the quarks that strike the off-shell boson (left understood), as well as on the heavy-quark masses, as indicated. For simplicity, in the above equation, we have set both renormalization and factorization scales equal to the hard scale, µ 2 R = µ 2 F = Q 2 , so α s = α s (Q 2 ). Finally, the index a denotes the type of structure function under consideration. In our case, we will mostly consider either F 2 or the longitudinal structure function F L . Note that when Q 2 ∼ m 2 Z , we also have a non-negligible contribution from the parity-odd structure function F 3 . This contribution is not logarithmically enhanced at small-x in the massless case [49]; this remains true for neutral current DIS with massive quarks, however, in charged current DIS with massive quarks a small-x logarithmic enhancement appears. To our knowledge, the resummation of these logarithms in F 3 was never considered so far.
In order to diagonalize the convolution in Eq. (2.1) we consider Mellin moments of the structure functions where, as often done in the context of small-x resummation, we have introduced a non-standard definition for the moments of the coefficient functions and of the parton densities: dz z N P ij (z, α s ), (2.4) where P ij are the usual Altarelli-Parisi splitting functions. In particular, in momentum space, the leading logarithmic (LL) behaviour at small-z to any order n > 0 in perturbation theory is , i.e. the enhancement is at most next-to-leading logarithmic (NLL). Note that some care has to be taken when considering the LO contribution C (0) a,i , which is an O(α 0 s ) constant in Mellin space, and hence formally LL.
In this paper we construct the NLL resummation of DIS structure functions at small-x. In order to achieve this goal, we resum the first two towers of logarithmic contributions to the splitting functions, while we consider only the first non-vanishing tower of logarithmic contributions to the partonic coefficient functions. Note that in a previous work [61] we called this NLL resummation in DIS coefficient functions just LL, underlining the fact that it is the leading non-vanishing logarithmic enhancement. We refer to the counting of this previous work as relative logarithmic counting (i.e. relative to the leading non-vanishing logarithmic enhancement), while the one adopted here as absolute counting (i.e. using the overall powers of α s and 1/N ).
Henceforth, unless explicitly stated, we are going to work in Mellin space and leave the dependence on the N variable, as well as the other variables, understood. As common in studies of DIS, we perform a flavour decomposition. For our purposes it is enough to separate the structure functions into a singlet and non-singlet component (2.5) The singlet structure functions contain both gluon and quark (singlet) contributions where

7)
n f denoting the number of active quark flavours. The non-singlet structure function instead reads, see e.g. [49,70], where e i is the electric charge of the quark q i , i.e. its coupling to the photon. 1 Furthermore, we can collect the terms proportional to the singlet PDF, obtaining where we have defined the so-called pure-singlet coefficient function, C PS a,q = C S a,q − e 2 C NS a,q , (2.10) being e 2 the average squared charge. In this study we consider resummed structure functions matched to their fixed-order counterparts. To this purpose, we find useful to introduce resummed contributions, defined as the all-order results minus its expansion to the fixed-order we are matching to, ∆ n C = C res − n k=1 α k s C res,(k) , (2.11) where we are going to typically consider matching to NLO and NNLO structure functions, i.e. n = 1, 2, although DIS structure functions are also known, in the massless case, to three loops [72][73][74]. Moreover, many of the formulae involving the resummed contribution that we derive in what follows hold regardless of the order in perturbation theory we are matching to. For these cases we adopt a simplified notation that omits the index n: ∆ n C → ∆C. Note also that we will sometimes write the full resummed expression C res as ∆ 0 C.

Factorization schemes in presence of massive quarks
In the context of collinear factorization, mass collinear singularities due to massless quarks must be factorized, such that perturbative coefficient functions are finite. This is the case for the up, down and strange quarks that we always consider massless. For massive quarks, i.e. charm, bottom and top, collinear singularities are regulated by the quark mass and they manifest themselves as logarithms of the ratio Q 2 /m 2 . Despite the fact that the factorization of these contributions is not necessary in order to obtain finite cross-sections, if Q 2 m 2 , these logarithms become large and their all-order resummation becomes desirable in order to obtain reliable perturbative predictions. This resummation is obtained by factorizing the mass logarithms, in the very same way as done for the massless quarks.
Whether mass collinear logarithms are factorized or not for a given massive quark is a choice of factorization scheme. A scheme where the collinear logarithms for the first n f lightest quarks are factorized is called a scheme with n f active flavours. In such a scheme, collinear logarithms are resummed for light quarks, while they are treated at fixed order for heavy quarks, thus defining the (relative) concept of light and heavy. Note that while obviously a heavy quark is massive, a light quark could be either massless, e.g. up, down, strange, or massive, e.g. charm and bottom. Thus, n f can be 3, 4, 5. Note that n f = 6 is not phenomenologically relevant, especially in the context of DIS, because the hard scale of the process is at most comparable, but never much bigger, than the top mass and so the top will be always treated as a heavy flavour.
In several cases, mostly in the contest of PDF fits where data span a large range in Q 2 , it is convenient to define so-called variable flavour number schemes (VFNS), where the number n f of active flavours varies as a function of Q, such that the collinear resummation for a given massive quark is turned on only for scales where it is needed. More specifically, a VFNS is a patch of factorization schemes with subsequent values of n f , which switches from a value n f to the next one (n f + 1) at a given "heavy quark threshold" µ h , typically chosen of the order of the heavy quark mass. The relation between a scheme with n f active quarks and a scheme where the mass logarithms of the (n f +1)-th flavour are resummed, i.e. a scheme with n f +1 active flavours, is at the core of the construction of a VFNS and provides the ingredients to resum the collinear logarithms of the (n f + 1)-th flavour. For a recent review see Ref. [75].
Let us consider DIS structure functions in a scheme with n f (massless or massive) active flavours. As we increase the hard scale Q, we reach energy scales which are significantly bigger than the mass of the (n f + 1)-th quark flavour. In this situation the n f -flavour scheme is no longer appropriate, as potentially large collinear logarithms log(Q 2 /m 2 ) are left unresummed. A more reliable framework is then provided by a factorization scheme in which n f +1 flavours are considered active, i.e. they all participate to parton evolution, having factorized, and hence resummed, their collinear behaviour [76][77][78][79][80][81][82][83][84][85][86][87][88]. The PDFs in the two schemes are related by matching conditions in the n f + 1 scheme evolve through standard DGLAP equations, i.e. as they all were PDFs of massless quarks. Note however that the heavy quark PDFs, being generated by the matching conditions Eq. (2.12), have a purely perturbative origin, and depend effectively on the heavy quark mass.
The structure functions Eq. (2.9) can be written in either scheme To all orders in α s , the choice of factorization scheme is immaterial and the two expressions are identical. Truncating the perturbative expansion of the coefficients to any finite order makes the two expressions different by higher order terms. Requiring equivalence of the two expressions order by order allows to relate the various coefficients and to find the matching functions K ij . The coefficient functions in the n f scheme are computed in standard collinear factorization with n f active quarks, with the heavy quark(s) only appearing in the final state or through loops. In the n f + 1 scheme the coefficient functions generally differ as the collinear logarithms due to the heavy quark are also factorized. Their expressions, which include mass dependence, can be determined by the equality of the first and second line of Eq. (2.13): (2.14) and similarly for the quark coefficient functions. In the equation above, we have defined and similarly for the gluon to light-quark matching function K qg , since these functions do not depend on the specific flavour but only on whether the final quark is light or heavy.
Eq. (2.14) and its quark counterparts can be solved to express the coefficient functions in the n f +1 scheme in terms of those in the n f scheme and the matching functions K ij , as we shall see in the next section. Furthermore, requiring that resummation of collinear logarithms due to the n f + 1 flavour is achieved in the n f + 1 scheme allows us to derive expressions for the matching functions as well. We will come back to this point in Sect. 2.1.2. For a detailed and general discussion, not limited to small x, see e.g. [75,88].

Heavy-flavour schemes at small-x
We now focus our discussion of heavy-quark factorization scheme to the contributions that are enhanced in the small-x regime. This topic has been discussed at length in Refs. [42,63], where explicit resummed results were presented in the DIS factorization scheme. Here, we extend the results to MS-like schemes.
As previously mentioned, the resummed contribution ∆C can be either seen as a contribution to the singlet or to the pure singlet. At the accuracy we are considering here, it always comes from a gluon ladder which ends with a quark pair production, one of which is struck by the photon, as depicted in Fig. 1. In the case of a quark initiated contribution, the quark immediately converts to a gluon, which then starts emitting. Therefore, the resummed contribution to the singlet coefficient functions, denoted ∆C a,i , always has the form (in the n f scheme) where ∆c a,i (m q k ) is the resummed contribution in the case of a light active flavour being struck by the photon, and ∆ c a,i (m q k ) is the resummed contribution in the case of a heavy flavour being struck by the photon. 2 Recall that being active does not necessarily imply being massless and indeed both massless and massive flavours contribute to ∆c a,i . Crucially though, in this contribution collinear logarithms are factorized and resummed and, consequently, the zero mass limit of ∆c a,i is finite. On the other hand, ∆ c a,i (m q k ) only contains massive quarks and no resummation of mass logarithms has been performed. Thus, the massless limit of this type of contributions is logarithmically divergent. In some simplified approaches, the mass of the heavy-quark is immediately neglected once it becomes active: this leads to what is sometimes called a zero-mass variable flavour number scheme (ZM-VFNS). Note that the massless contribution is identical for each massless quarks, so in a ZM-VFNS we would have n f k=1 e 2 k ∆c a,i (0) = e 2 n f ∆c a,i (0), i = g, q. (2.17) For this reason, a factor n f is usually included in the definition of the massless singlet coefficient function, see e.g. Refs. [49,61]. Here instead, we wish to retain the mass dependence of the active flavours, if present. To this purpose, we adopt a factorization scheme akin to S-ACOT [78,79] or FONLL [84] (which are formally identical [75,88]) in which the mass dependence is retained in the coefficient functions. We now perform a logarithmic counting on Eq. (2.14). None of the matching functions are LL, with the exception of the LO diagonal components, which are all equal to 1 at O(α 0 s ), K 2 In charged-current DIS, where the photon is replaced by a W boson, the quark flavour changes after hitting it, and so does its mass. Therefore, in this case, the coefficient functions ∆c a,i and ∆ c a,i would also depend on the mass of the outgoing quark. ), which start beyond NLL, as they are given by conversions of gluons or quarks into quarks with the participation of the heavy quark, i.e. suppressed by at least two genuine powers of α s . Note that, for simplicity, we are not indicating in ∆K ij the label [n f ] , but we emphasize its (logarithmic) dependence on the heavy quark mass.

∆C
Note that C NS,(0) a,q in Eq. (2.18) is the one in the n f + 1 scheme, so it is, in principle, an unknown of the problem. However, the requirement that in the n f + 1 scheme the resummation of collinear logarithms is achieved, forces it to be equal to the value computed in the limit where the heavy flavour is massless, up to possibly power-behaving mass corrections. This mass dependence can be arbitrarily fixed to be zero, as the original set of simultaneous equations, Eq. (2.14), is undetermined: there are two more unknowns than equations [75,88]. This choice is the one leading to S-ACOT/FONLL [78,79,84] and TR [80,81] schemes, where we have In alternative approaches, like ACOT [76,77,85] or, equivalently, a new incarnation of FONLL [75,89], denoted here as FONLL IC , which has been introduced to account for a possible intrinsic component of the charm PDF, the incoming heavy quark is treated as massive and therefore the mass-dependence is fully maintained in the coefficient function: . (2.20b) Note that the massless limit of these expressions reduces to the massless result, Eq. (2.19). We stress that for most applications the simpler S-ACOT/FONLL option is formally, and practically, as good as the more complicated ACOT/FONLL IC approach and hence we focus on it in the following. However, care must be taken when describing the charm structure functions in the case in which the charm PDF is fitted. In this case, the two approaches may lead to sizeable differences and the use of ACOT/FONLL IC might be advisable. We will come back to this point in Sect. 2.2.3. We now consider again the decomposition Eq. (2.16). In the n f + 1 scheme it simply becomes where ∆c a,i (m qn f +1 ) are the small-x resummed contributions to the production of a heavy quark (pair) in the scheme in which such heavy quark participates to the parton dynamics and its collinear logarithms have been factorized. We now plug Eqs. (2.16) and (2.21) into Eq. (2.18). All terms involving the lightest n f flavours and the heaviest n f + 2, . . . , 6 flavours cancel out, leaving only a relation between the coefficient functions for the n f + 1 flavour: The squared charge is the same in all terms and cancels out, thereby suggesting that this result will hold in general for more generic couplings, such as when the Z-boson exchange plays a role, as will be discussed in Sect. 2.1.3. We then find the (collinearly factorized) massive coefficients in the n f + 1 scheme for F 2 , assuming Eq. (2.19), (2.23) and for F L Given the matching function ∆K ij at NLL accuracy, the above results completely fix the relation of the NLL resummed contributions to the coefficient functions in n f and n f + 1 schemes. In particular, Eq. (2.24) shows that the resummed contribution to F L is the same in either schemes.

Computation of matching functions
In the derivation above the matching functions K (or ∆K ij ) are assumed to be given as an input to the computation of coefficient functions in the n f + 1 scheme. However, the very same derivation also allows us to construct the matching functions themselves. This is true in general (see e.g. Ref. [75]), but we focus on the small-x limit for simplicity.
The key observation is that, after scheme change, the coefficient functions in the n f + 1 scheme must not contain anymore the collinear logarithms associated to the (n f + 1)-th flavour. This is possible only if the matching functions subtract such collinear logarithms from the massive coefficients ∆ c a,i (m qn f +1 ), such that the massless limit m qn f +1 → 0 of the coefficient functions in the new scheme ∆c a,i (m qn f +1 ) is finite. If we further require that the n f and n f + 1 scheme are both of the same type, e.g. MS-like, we also need to impose that the massless limit of the massive coefficient ∆c a,i (m qn f +1 ) is just what we would have computed if the (n f + 1)-th flavour were massless: lim (2.25) (Note that, at the logarithmic accuracy we are interested in, ∆c a,i (0) are the same in the n f and in the n f + 1 schemes.
3 ) The massless limit Eq. (2.25) ensures that in the n f + 1 scheme the collinear logarithms are properly factorized into the PDFs and resummed through DGLAP. It also fixes the "constant" (i.e., non mass dependent) part of the function ∆c a,i (m qn f +1 ). It does not tell us anything about the power corrections in m qn f +1 /Q in the n f + 1 scheme, which have to be determined by the matching procedure. The results Eq. (2.24) show that, since F L has no collinear singularities at NLL, the massive coefficient in the n f + 1 scheme smoothly approaches the massless one at large Q, without any scheme change to be applied. On the other hand, in the case of F 2 , Eq. (2.23), which contains collinear singularities in the massless limit, the scheme-change effectively subtracts the matching function, and the difference will smoothly tend to the massless coefficient for Q m qn f +1 . We can exploit the last consideration to derive the desired resummed expressions for the matching functions ∆K hi , i = g, q. Indeed, Eq. (2.23) can be inverted to give ∆K hi in terms of the massive coefficient functions ∆ c 2,i (m), which are known, and of the coefficient functions ∆c 2,i (m) in the n f + 1 scheme, the massless limits of which, ∆c 2,i (0), are known. We can therefore consider the massless limit m qn f +1 → 0 of this expression, keeping in mind that ∆K hi and ∆ c 2,i (m) are separately logarithmically divergent, and therefore the massless limit has to be intended as setting to zero power-behaving contributions while keeping the logarithms finite. Assuming, or better choosing, that ∆K hi contain only logarithms of the mass m qn f +1 or mass independent terms, we then find where the limit has to be intended as a formal expression as described above. Note that including power-behaved mass dependent terms is possible and would simply change the form of ∆c 2,i (m) through Eq. (2.23), as well as the PDFs in the n f + 1 scheme according to Eq. (2.12). However, this "factorization" of power behaved contributions is beyond the control of the collinear factorization framework, so the simplest choice Eq. (2.26) is perfectly acceptable and therefore universally adopted. The definition Eq. (2.26) also shows that matching functions at NLL satisfy colour-charge relations. Indeed, it is known that both massless and massive coefficient functions in the quark and gluon channels are related at NLL by 4 for the coefficient functions in the n f + 1 scheme.

Generalization to neutral and charged currents
The previous results have been derived in the case of photon-exchange DIS. We now comment on their generalization to the more general case of neutral current, where also the Z boson contributes, and the case of charged current, where the exchanged boson is a W . Including a Z boson in the discussion is rather trivial. What changes is just the coupling to the quark, which is different to that of the photon both in the Z-only exchange and in the photon-Z interference. In Z-only exchange both vector and axial couplings contribute, which we denote generically as g V and g A . If we consider only massless quarks, the sum of the squares of each coupling, g 2 V + g 2 A , factorizes. The only subtlety regarding Z exchange is that in the massive case there is a contribution which is not proportional to the sum of the squares of the vector and axial couplings. However, in this case, one can still factor out g 2 V + g 2 A , but a contribution proportional to g 2 A /(g 2 V +g 2 A ) appears, which is not present in the photon (or photon-Z interference) case. This term is not problematic, as it vanishes in the limit of massless quarks. Therefore, the whole construction of the previous sections, and in particular Sect. 2.1.1, remains unchanged, provided the photon couplings e 2 k are replaced with the more general coupling (g 2 V + g 2 A ) k for each quark q k . The charged-current case is less trivial. In this case, the quark flavour changes after interacting with the W and this results in two additional complications. Firstly, the quark mass changes. In fact, for all practical applications, one of the quarks interacting with the W can be considered as massless. Indeed, ignoring the top quark which is too heavy to give a contribution at typical DIS energies, the only two massive quarks are the charm and the bottom, but their interaction is suppressed by the V cb CKM entry and can thus be neglected. Hence, the only remaining combinations involve either a massless and a massive quarks, or two massless quarks. The second complication arises due to the fact that the final state, composed by a quark q and an anti-quarkq , is not self conjugate (as it is in the neutral-current case, where the final state contains qq). This means that the non-singlet coefficient in Eq. (2.13) is different for q andq, depending on the process. In particular, this implies that when using Eq. (2.14) to derive Eq. (2.18), only either the heavy quark q n f +1 or the heavy anti-quarkq n f +1 contributes: thus, the collinear subtraction, where present, contains only half of the matching function. Additionally, for the computation of the parity-violating F 3 structure function (which contains a singlet contribution in the massive charged-current case), the non-singlet LO coefficient has opposite sign depending on whether the initial state parton is a quark or an anti-quark, specifically C where we are denoting with Q the heavy massive quark and with q the companion massless quark appearing in the final state. The massless limit of the latter is just zero, ∆c CC 3,i (0) = 0, since in the massless limit F 3 is non-singlet and therefore it does not contain logarithmic enhancement at small x (see App. A for further details).

Small-x resummation of coefficient functions and matching functions
We are now ready to discuss the actual resummed expressions for coefficient functions both in the massless and massive case. Additionally, using Eq. (2.26), we also determine the all-order behaviour of the matching functions. The all-order behaviour of partonic coefficient functions is obtained using the k t -factorization theorem. In this framework one computes the gluon-initiated contribution to the process of interest, keeping the incoming gluon off its mass-shell. To the best of our knowledge, only the off-shell coefficient functions that are necessary to perform the resummation of the photon-induced DIS structure functions have been presented in the literature, both in the case of massless and massive quarks, see e.g. [44,45,90]. However, in this study we want to resum all DIS coefficient functions, both the neutral-current and charged-current contributions. Therefore, we perform a general calculations of DIS off-shell coefficient functions considering the coupling to the electro-weak bosons and, where relevant, the interference with the photon-induced contributions. The calculation is detailed in Appendix A, where we collect all the relevant results for the off-shell DIS coefficient functions. These results have been implemented in the public code HELL, version 2.0, and are also accessible through the public code APFEL [91], to which HELL has been interfaced, which directly constructs resummed structure functions.

Resummed coefficient functions
The resummation of the massless coefficient functions was originally performed to NLL in Ref. [49]. Following Ref. [38], in our recent work [61] we have also included formally subleading, but important, contributions such as the ones related to the running of the strong coupling. In our approach, the partonic massless coefficient function C a (N, ξ, 0, 0), which is calculated with an incoming off-shell gluon (see App. A for details about our notation), is convoluted with an evolution operator U (N, ξ), and γ + is the small-x resummed anomalous dimension. Note that the above expressions hold in a factorization scheme denoted Q 0 MS [47,49,92,93], which differs from MS at relative O(α 3 s ). In the context of small-x resummation this scheme is preferred because it leads to more stable results [38]. Furthermore, since the explicit N dependence of the off-shell partonic coefficient function is subleading, we find advantageous to work at NLL with its N = 0 moment. In the resummed expression of the F 2 contribution, the subtraction term U qg appears. Its role is to cancel the collinear singularity of C 2 (0, ξ, 0, 0). Its expression reads [61] where γ qg is the resummed qg anomalous dimension. All ξ integrals extend to ∞ and start from the position of the Landau pole, ξ 0 = exp −1 αsβ0 . In ξ = ξ 0 the evolution function U is supposed to vanish (this was e.g. a condition for neglecting a boundary term when integrating by parts in Ref. [61]). However, due to subleading contributions, this is not always true in our practical construction. While the induced effect is subleading, this fact is undesirable: we discuss in Appendix B.1 how we now deal with this issue.
We can apply the same procedure to the case of a heavy (non-active) flavour. We start considering neutral currents. Note that in this case the mass of the quark acts as a regulator and no subtraction term U qg appears, where we have defined ξ m = m 2 /Q 2 (see again App. A for the precise definition of the off-shell coefficient and its arguments). The massive coefficient functions entering the above formula have been computed a long time ago [45,94] (see also Ref. [95]). We report them in Appendix A, where we have also recomputed them in a more general set-up which covers the full neutral-current case in which the mediator can be a Z boson thus finding a new contribution proportional to the axial coupling, and the charged-current scenario in which the mass of the quark changes after interacting with the W boson. We have already commented in Sect. 2.1.3 that for physically relevant processes only one of the quarks involved in charged current DIS is massive, and the other can be treated as massless. In this case the resummed coefficients read Here, since one of the two quarks involved is massless, we need massless collinear subtractions, implemented through U qg , to take care of the collinear singularity. Since only one out of two diagrams contains the singularity, there is a factor 1/2 for each subtraction. Each subtraction further multiplies the LO non-singlet diagram evaluated in N = 0, corresponding to the process q + W → Q (or its conjugate) with massive Q, which has non-trivial mass dependence for F L and F 3 (and non-trivial sign for F 3 ). For F L in particular, this term vanishes in the massless limit, consistently with Eq. (2.31b). In the last equation we are treating separately the cases in which the final state contains a massless quark plus a massive anti-quark and its conjugate process, for the same reason discussed in Sect. 2.1.3. Note that, according to the notation defined in App. A where the third argument of the off-shell coefficient is the mass (squared divided by Q 2 ) of the anti-quark and the fourth the mass of the quark in the final state, the arguments are swapped in the two cases. Effectively, for F 3 swapping the arguments changes the sign, so the difference between the two cases is just an overall sign. For F 2 and F L , instead, the coefficients are symmetric for final-state charge conjugation, and therefore the result does not change when swapping the arguments. We stress that in the massive case the partonic coefficients include non-trivial theta functions which restrict the available phase space. This is originally encoded in the N dependence of the offshell coefficients, which we loose when setting N = 0. As these theta functions are very physical, it is important to restore them. Details on how this is implemented in our resummed results are given in Appendix B.2.
The massless ξ m → 0 limit of all the off-shell coefficients is finite, because the off-shellness ξ regulates the collinear region, and gives the massless off-shell coefficients entering Eqs. (2.31) (and zero for F 3 ). However, while the ξ m → 0 limit for F L gives automatically the massless result, as it must according to Eqs. (2.24), (2.30b), for F 2 and F 3 one further needs to subtract the matching condition, Eqs. (2.23), (2.30a), (2.30c). Accordingly, the resummed massive coefficient functions for the massive active flavour are given in neutral current by and in charged current by (again, the last equation is split in two depending on whether the final state is qQ or Qq, the difference being an overall sign). Comparison of Eq. (2.36a) (and equivalently Eq. (2.37a)) with Eq. (2.31a) shows that the massless limit does not commute with the on-shell limit in presence of collinear singularities. In particular, the ξ m → 0 limit applied to ∆ 0 c (CC) a,g (m), a = 2, 3 does not commute with the ξ integration.

Resummed matching functions
We can now use Eq. (2.26) to determine the resummed matching function ∆ 0 K hg (m). Using Eqs. (2.34) and (2.31a) we can write 6 where the last two terms are basically the commutator of ξ integration and massless limit. Computing this commutator in the general case is highly non-trivial. Therefore, we consider first the limit in which the coupling is kept fixed. In this limit the convolution over ξ becomes a Mellin transformation with moment M = γ s αs N , which is the the LL anomalous dimension, dual of the LO BFKL kernel. This Mellin transform can be performed analytically, so the massless ξ m → 0 limit can be safely taken afterwards. We have (see Appendix A.3) A non-trivial check of the above expression can be performed by considering its perturbative expansion, where h i are the (known) coefficients of the perturbative expansion of γ s U qg f.c.
= γ qg in powers of γ s [49,61]. The first coefficients read 14 9 . Note that the collinear pole 1/γ s cancels out in the sum. The second equality is then obtained replacing γ s = αsC A πN + O(α 2 s ). The above expression can be compared with the fixed-order results presented in Ref. [82]. The α s term corresponds to the Mellin transform of the NLO result computed in N = 0, and the α 2 s term correctly reproduces the leading singularity of the NNLO contribution. Checking the above result one order higher in perturbation theory is less trivial because at this order we start to become sensitive to the choice of factorization scheme. After taking into account the conversion from Q 0 MS to MS, which affects the second term in Eq. (2.39), we find full agreement with the high-energy limit of the α 3 s result [96,97]. We can now restore the running-coupling effects in the resummation from the fixed-coupling result by computing where K hg is obtained as the inverse Mellin transform of its fixed-coupling counterpart, second term in Eq. (2.39). The computation of this inverse Mellin transform is done in App. A.3, and its derivative is given by Clearly, by construction, Eq. (2.41) with Eq. (2.42) reproduces the correct result in the fixedcoupling limit, Eq. (2.39). Also, it clearly includes the correct resummation of the subleading running coupling effects, as the form Eq. (2.41) is the standard expression for such resummation [61]. Eq. (2.42) is a new result. It allows to resum the matching conditions in MS-like schemes with full inclusion of running-coupling effects.

Matching to fixed order and construction of the VFNS: FONLL as an example
We conclude the section by giving some details on how the resummed results presented above are combined with the fixed-order contributions to construct a VFNS. First, we have to subtract from the resummed coefficient functions their expansion up to the order we want to match onto. The O(α s ) and O(α 2 s ) contributions to the resummed matching function are explicitly given in Eq. (2.40) in the fixed-coupling limit, but up to this order they are identical to their running coupling counterparts, and so they can be used straight away to construct ∆ 1 K hg (m) and ∆ 2 K hg (m). The same equation shows in the first line the expansion of The actual construction of a VFNS is more delicate. Indeed, there are at least two degrees of freedom that have been exploited in the literature to construct different incarnation of VFNSs (at fixed order). One degree of freedom is related to the inclusion of undetermined (by the matching conditions) power-behaving mass dependent contributions in some coefficient functions, as already discussed in Sect. 2.1.1. The second degree of freedom is related to how to combine the various ingredients at a given finite perturbative order. The approach adopted so far in our construction can be identified with a plain (i.e., without any χ-rescaling [98]) S-ACOT construction, with a canonical perturbative counting based on explicit powers of α s (at fixed α s /N for resummed contributions). This is equivalent to a plain (i.e., without damping [84]) FONLL construction, even though in FONLL the various ingredients are combined together with a different philosophy. In the following we will briefly review the FONLL construction, which is implemented in the APFEL+HELL package, with which our results can be directly used for resummed DIS phenomenology.
The FONLL approach [83] is a standard combination of fixed-order and resummed contributions, in which these two ingredients are simply summed up and the double counting between them subtracted. In the FONLL case, the distinction between "fixed order (FO)" and "resummation (NLL)" refers to collinear logarithms due to massive quarks. In DIS, the FONLL construction [84] of structure functions, assuming a single heavy quark q n f +1 with mass m, is performed as (m) is the fixed-order (called massive) contribution, in which the collinear logarithms are not resummed and which retains the full mass dependence of the heavy quark, F (0) is the resummed (called massless) contribution, computed assuming that the heavy quark is massless, and thus where the (singular) collinear logarithms are resummed, and finally F d.c. a (m) is the double-counting term (called massive-zero), which can be either seen as the fixed-order expansion of F [n f +1] a (0) or the "massless limit" (in which divergent terms are kept finite) of F a (m) can be seen as the resummed contribution to be added to the fixed order to resum the collinear logarithms, or equivalently F a (m) can be interpreted as the power-behaving mass corrections to the (resummed) massless calculation.
According to our nomenclature (extended to the structure functions) the FONLL result is just the massive result in the n f + 1 scheme, i.e.
Thus, the double counting term, which is the only non-trivial ingredient in Eq. (2.43), is given by for NC, and similarly for CC (with an extra factor 1/2 multiplying ∆K hi ). The discussion so far does not add anything to the results presented in the previous sections. However, having now the small-x resummation for each individual ingredient appearing in Eq. (2.43), we can also consider the version of FONLL which includes a damping on the resummed contribution, such that the resummation smoothly turns off at the scale Q = m. This variant is often used in PDF fits, and effectively corresponds to damping the collinear subtraction term −C NS,(0) a,q ∆K hi (m) in our resummed coefficients ∆c a,i (m).
A word of caution is needed when discussing the perturbative counting. The canonical counting would consist in including all contributions at O(α s ) for NLO and all contributions at O(α 2 s ) at NNLO (and so on). However, a non-standard counting is usually adopted at NLO (e.g. in NLO NNPDF fits), where the massless contribution F  [84]. When this particular perturbative counting is adopted, the double-counting piece must be computed with care. In particular, only the definition of F d.c. a (m) as the fixed-order expansion of F s ) gives the correct result. As far as small-x resummation is concerned, one has to use ∆ 2 c a,i (m) for the massive part and ∆ 1 c a,i (0) for the massless part, while for the matching at the heavy quark threshold in DGLAP evolution ∆ 1 K hi (m) is to be used. For the double-counting part, being it the expansion of the massless, the use of ∆ 1 c a,i (0) and ∆ 1 K hi (m) in Eq. (2.46) is needed. It can be explicitly verified that in this way the "resummed contribution" F s ) and subleading at small x. Finally, we discuss a variant of the VFNS where the mass of the heavy quark is retained in all coefficient functions in the n f + 1 scheme. This is the original ACOT [76,77,85] construction, and it also corresponds to the variant FONLL IC proposed in Refs. [75,89] to account for a possible intrinsic component of the charm PDF. Following the latter references, we define (2.48) to be the term needed to upgrade S-ACOT/FONLL to ACOT/FONLL IC (ignoring damping, rescaling, etc.). The small-x resummation of this term can be simply obtained by computing the difference between the resummation obtained with massive NS coefficients, Eq. (2.20), and the one obtained with massless NS coefficients, Eq. (2.19). We thus have (we apologize with the Reader for the awkward notation) which are the resummed contributions to the singlet coefficient functions δc a,i making up δF a for neutral current and charged current respectively. Note that the massive NS coefficients, which have a non-trivial dependence on N , can be computed in N = 0 in Eqs. (2.49), as the N dependence is a subleading effect as small x.

A new approach to running-coupling resummation in DGLAP evolution
In the original ABF construction [33][34][35][36][37][38], which we followed with minor modifications in our previous work [61], the resummation of the anomalous dimension γ + (the largest eigenvalue of the singlet sector) is performed through the exploitation of the duality relation between DGLAP and BFKL evolution kernels, improved with symmetrization of the latter and the imposition of exact momentum conservation. This result is usually referred to as double-leading (DL) resummation. However, it was realized long ago [35] that running coupling corrections to fixed-order duality give rise to subleading terms which potentially spoil the perturbative stability of the result. Therefore, despite their formally subleading nature, the resummation of these effects is of utmost importance in order to obtain stable and reliable resummed anomalous dimensions. Additionally, the resummation of these terms changes the nature of the all-order small-N singularity, converting a square-root branch-cut into a simple pole. Therefore, the resummation of these contributions, known as running-coupling (RC) resummation, is usually added to the DL result.
The RC resummation can be obtained by solving the BFKL equation with full running coupling dependence (see e.g. [39][40][41]), and then deriving from the solution (which is an eigenvector PDF) its anomalous dimension. If we were able to perform this procedure with the full DL BFKL kernel, the resulting anomalous dimension would just be the final result. However, solving such equation analytically is not possible, due to the complicated all-order α s -dependence and the non-trivial M -dependence of the DL kernel. In some approaches, e.g. [28][29][30][31][32], the equation is thus solved numerically, and the resummed anomalous dimension derived in a numerical way. Instead, in Refs. [37,61] an approximate analytic solution, in which both α s -and M -dependencies of the kernel are simplified, was constructed and added to the DL anomalous dimension, subtracting the appropriate double counting.
We find this second approach more convenient, and we keep adopting it here. However, in this section we critically review the approximations used in Refs. [37,61] and propose a new way of constructing and approximating the kernel from which the RC solution is computed. Our approach has various advantages, from purely practical ones related to the numerical implementation to most serious ones related to the physical nature of the solution and its α s dependence.

The choice of the kernel
The core of small-x resummation of the largest eigenvalue γ + is encoded in the duality between the DGLAP and BFKL equations. Imposing that the corresponding eigenvector PDF is solution to both equations requires the duality relation where χ DL (M, α s ) is the BFKL kernel and γ DL (N, α s ) the DGLAP anomalous dimension. The knowledge of the BFKL kernel at N k LO provides by duality all the N k LL contributions in the anomalous dimension, and vice versa. The name DL comes from the fact that both kernels are supposed to contain their fixed-order part at N k LO, thus implying that they also contain (by duality with each other) all N k LL contributions. Therefore, dual DL kernels obtained with fixed N k LO in both, usually denoted DL-N k LO, are matched results of the form N k LO+N k LL, and so they are both double (next-tok ) leading order and log. The actual DL kernel and anomalous dimension further contain additional ingredients (from symmetrization and momentum conservation) which are required to make the result perturbatively stable. The duality relation Eq. (3.1) assumes that the strong coupling α s does not run, namely it is Q-independent. When the running of α s is taken into account, the duality equation receives additional corrections. If these corrections are included perturbatively, new singularities appear which make the result perturbatively unstable. For instance, at NLO+NLL one should include a purely NLL term of the form of a LL function of α s /N times an overall factor α s [35] The new singularity is obtained when χ 0 (M ) in the denominator vanishes, i.e. in M = 1/2. Higher-order corrections will have larger powers of χ 0 (M ) in the denominator, leading to a perturbatively unstable singularity for N such that γ s (α s /N ) = 1/2, i.e. N = α s χ 0 (1/2). This instability goes away if RC corrections are included to all orders in the running-coupling parameter β 0 (i.e., at NLO+NLL one should include a term of the form of a LL function of α s /N times a function of α s β 0 to all orders), and the various singularities sum up to a simple pole, whose position is perturbatively stable. This is the main motivation for including the RC corrections to all orders. RC corrections are resummed by solving the BFKL equation with the DL kernel χ DL in which α s is not fixed but is running. When the coupling runs, α s (Q 2 ) becomes a differential operator , Q 0 being a fixed scale, and where we have assumed 1-loop running. In principle, there can be α s computed at different scales in the kernel, but one can always rewrite it as α s (Q 2 ) by evolving to that scale and expanding the relevant evolution factors. In this way, theα s operators are placed on the left of all the M -dependent terms of the kernel, and act on everything to the right. This ordering is the one chosen by ABF to derive their solution of the RC equation. Two observations are now in order.
• While at DL-LO all the running coupling evolution factors are higher orders, and α s can be set equal to α s (Q 2 ) in all terms without modifying explicitly the kernel, at DL-NLO changing the argument of each α s to α s (Q 2 ) in all LO contributions produces terms which are formally NLO and have to be included in the kernel. Because of the symmetry properties of the BFKL ladder, the DL-NLO kernel (see Eq. (4.14) later) does not correspond, by construction, to a kernel in which all α s 's are α s (Q 2 ). Therefore, at NLO, the form of the kernel in which all α s 's are α s (Q 2 ) differs from that of the DL kernel by NLO terms. For this reason, a different kernel was used for DL and RC resummation at NLO in Refs. [37,38,61].
• Once all powers of α s are computed at Q 2 , which is equivalent to say that all powers ofα s have been commuted to the left, the running coupling equation cannot be solved directly, because it is in principle a differential equation of infinite order. Therefore, in ABF a linear approximation of the kernel, in which at maximum one power ofα s is retained and all others are frozen to α s (Q 2 0 ), is used. The fact that the complicated all-orderα s dependence of the kernel is approximated with a linear one may seem too crude. However, the goal of the RC resummation is to resum to all orders a class of terms, behaving as powers of α s β 0 times a LL function of α s /N , which originates from 1-loop running at lowest order. The NLO contribution to the BFKL kernel would produce corrections which are of order α s (α s β 0 ) n (α s /N ) k for all n, k, and therefore beyond the formal accuracy we aim to. This shows that the linear approximation suggested by ABF is in fact sufficient to the scope of this resummation.
In fact, this argument also suggests that using a NLO kernel for RC resummation which differs from the DL one is unnecessary. Indeed, the ingredients which determine the leading RC corrections to all orders are contained in the LO part of the kernel. Thus, correcting the kernel by NLO terms will change subleading RC contributions which we do not aim to resum. Therefore, for the accuracy we are interested in, the NLO part of the RC kernel is immaterial, and there is no reason for using a different kernel for RC resummation and DL resummation.
Using the same kernel for DL and RC resummations, as we now suggest, has an important consequence: we do not need any more the infamous function γ match introduced in Ref. [37] and used later in Ref. [38,61] to cure the mismatch of singularities between the DL part and the RC part of the result. This function was needed to effectively remove the square-root branch-cut of the DL solution, since the subtraction of the double counting term between DL and RC resummations does not cancel it, exactly because two different kernels were used. What we have realized only with this work is that the function γ match is not subleading as was originally claimed [37], but rather it is NLL and therefore ruins the formal accuracy of the NLO+NLL result. We give more details about this function in Sect. 3.4. Having understood that the linearα s approximation is a good approximation, and that as such the same BFKL kernel can be used for both DL and RC resummations, in order to be able to solve the RC BFKL equation we further need to specify the M functional form. The approximation adopted by ABF, which we followed in Ref. [61], is a quadratic approximation around the minimum of the kernel. Indeed, the minimum encodes, by duality, the information on the leading singularity, and it is therefore sufficient to accurately describe the kernel and to perform the RC resummation. However, we argue that this quadratic approximation has subtle undesired properties which makes it not ideal for our purposes. For instance, the α s -expansion of the resulting anomalous dimension contains half-integer powers of α s . This is a direct consequence of the fact that a polynomial kernel, such as this quadratic approximation, is non-physical. Thus, here we propose a different approximation, which is physically motivated and which leads to an expansion in integer powers of α s . We discuss this new approximation, denoted collinear approximation, in the next section.

Solution of the RC differential equation in the collinear approximation
We are now going to derive the solution of the RC BFKL equation in the linearα s approximation and collinear M approximation. The starting point is the on-shell BFKL kernel in symmetric variables [37] that we denote simply χ(M, α s ), whoseα s dependence is approximated as where prime denotes derivative with respect to α s . It is important to observe that this approximation includes an O(α 0 s ) term,χ, which is not physical and not present in the original kernel which is of lowest order O(α s ). For this reason the kernel Eq. (3.3) does not go to zero asα s → 0. One could therefore consider another linear approximation, which does go to zero asα s → 0, but does not reproduce the exact derivative inα s = α s . In fact, both approximations are equally valid for our purposes, as they would be identical (and exact) in the case of a LO kernel χ(M,α s ) =α s χ 0 (M ). In the following, we will use the first approximation, Eq. (3.3), to derive our results, which can then be easily translated into the second approximation, Eq. (3.4), by simply lettingχ = 0, χ = χ/α s . We stress that this simple translation rule could not be applied in the original ABF solution with quadratic kernel, as the limitχ → 0 is not trivial in that case. The (homogeneous) RC BFKL equation with kernel Eq. (3.3), from which the resummed anomalous dimension can be derived, is given by [37] where f (N, M ) is the double Mellin transform of the eigenvector PDF. Assuming 1-loop running, and taking the logarithmic derivative of the solution, we arrive at the anomalous dimension where M 0 and M 1 are free parameters which must be in the physical region 0 < M < 1 (they can be conveniently chosen to be equal to each other, and equal to the position of the minimum). In order to compute the integrals analytically, we need to specify the form of the kernelsχ and χ . In the ABF construction, a quadratic approximation around the minimum of the actual kernel was considered, where M min differs in general from 1/2 by terms of O(α s ). The polynomial form of the quadratic kernel is non-physical, as for instance the inverse Mellin transform of the n-th power of M corresponds to the n-th derivative of a δ function of k t in momentum space. A better approximation, which we propose here, is inspired by a collinear plus anti-collinear approximation of the kernel [93,99] generalized to account for a minimum which is not in 1/2: Expanding around its minimum M = M min we find which leads to the identifications such that the collinear kernel incorporates exactly the same information as the quadratic kernel. Therefore, from the point of view of the accuracy of the approximation, the new collinear kernel is as good as the old quadratic one. However, as its form resembles the leading collinear and anticollinear poles of the actual kernel, it performs better than the quadratic one, and leads to a solution with better features, as we shall now see.
To compute the solution of the BFKL equation we need to specify theα s dependence of the collinear kernel Eq. (3.8). For A(α s ) and B(α s ), we use the very same linear decomposition Eq. (3.3), while the position of the minimum M min (α s ) will not be considered as an operator. 7 The integrals in Eq. (3.6) can be computed easily by noticing that the functional form of the integrand is identical to that of the quadratic kernel used by ABF, for which the solution is already known. Specifically, for quadratic kernel, the kernel-dependent part of the integrand of Eq. (3.6) is given by while for collinear kernel we find having used in the last step Eq. (3.10). By direct comparison, we find the translation rules Hence, the final solution is given by [37,61] We immediately observe that the limitc,κ → 0 of these expressions is finite and trivial, so the solution in the approximation Eq. (3.4) is a trivial limit of this solution. This is in contrast with the analogous solution with a quadratic kernel, whoseχ → 0 limit is not trivial and leads to a solution in terms of Airy functions. This represents a first advantage of using the collinear kernel with respect to the quadratic one. Eq. (3.14) with Eq. (3.15) represents a new result with respect to the old "Bateman" RC solution of Refs. [37,38,61] and it generalizes it to the case in which the minimum is not in 1/2. In order to study the properties of this result, we start considering the saddle point expansion of Eq. (3.6), which is equivalent to an expansion in powers of β 0 of Eq. (3.14). This expansion is also needed to identify the proper double counting with the DL part. We find This result shows a number of interesting features, especially when compared with the analogous expansion in the case of the quadratic kernel [37]. First, we note that at large N all terms go to zero except for the running coupling correction −β 0 α s , which is finite. This is in agreement with the fact that the starting kernel had a pole in M = 0, which by fixed-coupling duality leads to an anomalous dimension that goes to zero at large N . Note that this is in contrast with the case of the quadratic kernel, which diverges at large N as − 2N/κ, due to the absence in the first square root of the term +(N − c)/M 2 min in the denominator, in agreement with it being derived from a BFKL kernel quadratic in M . This term is crucial for another reason, as it makes the denominator of O(α 0 s ), while it would be of O(α s ) if the denominator were just κ/2, as it happens with the quadratic kernel. Having a denominator of O(α s ) produces an α s expansion of this result which contains half-integer powers of α s . Instead, the α s expansion of Eq. (3.16), and hence of Eq. (3.14), is perfectly acceptable with only integer powers of α s . These two differences represent two additional important benefits of using the collinear kernel rather than the quadratic one.

Construction of matched results
We recall that the solution Eq. (3.14), having been derived with an approximate M dependence, cannot be regarded as the full solution. Rather, it represents the all-order resummation of the β 0 terms which must be added to the DL result, after subtracting those contributions which are already included (and not approximated) in the DL result. In the following we will thus focus on the combination of our RC resummed anomalous dimension with the DL one, also providing the α s expansions of the RC contributions which will be needed in Sect. 4 for matching (N)LL resummation to (N)NLO.
When matching the RC resummation to the DL-LO result, the first three terms of the singular expansion Eq. (3.16) have to be subtracted (the first two because they are LL, and the third because it is of O(α s ), and hence already included in the DL-LO). After subtraction we have the expansion where κ 0 and c 0 are the derivatives of κ and c computed in α s = 0, and are given by The careful Reader might wonder what are the consequences of starting with a BFKL kernel in symmetric variables. Indeed, when combined with the DL result, the RC result must be translated to DIS (asymmetric) variables. This amounts to adding N/2 to the RC anomalous dimension. However, such a term is automatically subtracted in the construction of the RC contribution to the DL result, ∆ DL-LO γ rc (N, α s ) [100]. Therefore, the latter object is insensitive to the change of variables.
For NLL resummation the RC result must be matched to DL-NLO, so we further need to subtract the O(α 2 s ) of γ rc . However, we observe that at O(β 0 ) the expansion of γ rc , Eq. (3.16), contains terms which are formally NLL, and specifically given by α s β 0 times a LL function. These terms should be already included in the DL-NLO result. In fact, contributions of this form originate from running coupling corrections to the duality relation [35], Eq. (3.2), and are not automatically generated by fixed-coupling duality in the DL-NLO result. Rather, they have to be supplied to the DL-NLO result as an additive correction [37,61], where α s χ 0 (M ) is the LO BFKL kernel and γ s (α s /N ) its dual. The −1 in square brackets represents the subtraction of the double counting with the fixed-order part of the DL-NLO; after subtraction, this function starts at O(α 4 s ). Since the kernel used in the RC resummation is only approximate, the function γ rc does not correctly predict all the NLL contributions of Eq. (3.19). Therefore, Eq. (3.19) must be still added to the DL-NLO result, and the O(β 0 ) part of γ rc has to be considered as a double counting term with respect to ∆γ rc ss , and hence subtracted. Thus, for RC resummation matched to DL-NLO, we further need to subtract the fourth term of Eq. (3.16), where we have also written its α s expansion in the last line. We observe that, formally, only the O(α s β 0 ) times a LL function is really doubly counted, so in principle one could expand the last line of Eq. (3.20) at NLL, and remove the NNLL terms. However, we think that it is most consistent to also remove these spurious higher logarithmic order contributions. Interestingly, the expansions of both Eq.

Singularity mismatch
The At LO+LL, the kernel is the one obtained putting on-shell Eq. (4.2), so the singularity automatically cancels. 9 In Ref. [61] an intermediate result, denoted LO+LL , was introduced to perform the resummation of quark entries of the anomalous dimension matrix and of coefficient functions. This anomalous dimension is formally LO+LL, but uses the RC parameters of the NLO kernel such that the position of the leading singularity is the same as that of the NLO+NLL result. For this result, the cancellation of the branch-cut cannot take place. In order to cure the mismatch in the singularities of the DL-LO result and the RC result with NLO parameters we need a matching function γ match . Its form must be where the function γ m must reproduce the singular behaviour of the RC and DL parts, respectively, and the parameters c (N)LO , κ (N)LO and M (N)LO min are those obtained from the (N)LO kernel. For the case of collinear kernel, γ m may be simply given by the first two terms of Eq. (3.16). However, we have some latitude with the definition of the matching function as far as subleading corrections are concerned. We can exploit this freedom to define a matching which numerically has a very moderate effect. We find that the choice  [37,61] the kernel for RC resummation was constructed with a different α s ordering with respect to the DL one, which had different minima and thus created a singularity mismatch between the DL and RC anomalous dimensions, making it necessary the introduction of a matching function to cancel these singularities. We have already argued in Sect. 3.1 that using different kernels was not necessary, and we can actually use the same DL kernel also for RC resummation, thereby ensuring automatic cancellation of the square-root branch-cut. Therefore, in this work we no longer need to patch the NLO+NLL result with a matching function.
It is important to stress that, had we used two kernels for the DL and RC parts of the NLO+NLL resummation differing by O(α 2 s ) terms, the analogous matching function would have necessarily been NLL (as in LO+LL case), thus contaminating the result which could not be claimed to be NLL anymore. This is indeed the case for the one used in Refs. [37,61]. We have verified that it is not possible to modify the function γ m , Eq. (3.22), to make the function γ match NNLL without introducing new (uncanceled) singularities. To prove this statement, we consider a generalization of Eq.  (3.25) where the index in the derivatives indicates with respect to which argument the derivative is computed, and m 1 is the O(α s ) contribution to M min . The expansion of the square-root term at NLL depends on the NLO coefficients c 1 , κ 1 and m 1 . If the two kernel used for RC and DL resummations differ at O(α 2 s ), then these coefficients differ, and when building up the function γ match these NLL terms do not cancel among the two γ m 's. Thus, to only way to make the function γ match a purely NNLL object is to choose the functions η such that the NLL expansion of γ m vanishes. However, the expansion of the square-root term is singular at NLL in N = α s c 0 , so it is clear that in order to make γ match vanishing at NLL the derivatives of the function η, and thus the function itself, must be singular in N = c. But if this is the case, then the function γ m will contain additional singularities with respect to those which it is suppose to cancel. This violates our assumptions. We must thus conclude that it is not possible to cancel a NLO singularity mismatch and at the same time preserve NLL accuracy. This conclusion remains true also for the matching function used with the quadratic kernel. Therefore, the only way not to spoil the formal NLL accuracy of the NLO+NLL result is to use exactly the same kernel for the DL and the RC parts of the result: this is the main motivation for this choice, adopted in this work for the first time.
In the NLO+NLL result, there is a different singularity mismatch coming from the function ∆γ rc ss , Eq. (3.19), and the ∆ DL-NLO γ rc , Eq. (3.20). The latter exhibits explicitly a pole in N = c, which is different from an analogous singularity in the former, which is in N = α s c 0 , as one can easily verify from the definition. The singularities would be identical if the parameters of the RC kernel were those of the LO BFKL kernel, and would cancel in the sum. However, due to the higher orders contained in the parameters used to construct the RC kernel, the position of the singularity is shifted and the cancellation does no longer take place. We can solve the problem by introducing a new matching function to be added to the final result, which effectively replaces the singularity of Eq. (3.20) with Eq. (3.26). Being the singular contribution a NLL term, this matching function is formally NNLL, and therefore acceptable. However, as in the LO+LL case, it is convenient to subtract additional higher orders, such that the effect of this function is as moderate as possible. Our choice is where the last term (which is formally NNLL) ensures cancellation of a number of subleading contributions from the difference between the first two terms which could potentially spoil the accuracy of the result. Additionally, because of the last term, the function γ ss match starts at O(α 4 s ) and therefore it does not contribute to the α s -expansion of the NLO+NLL result to O(α 3 s ). Note that this singularity mismatch was present also in the original works using a quadratic kernel [37,61]; the problem there was solved by replacing by hand the singularity in ∆ DL-NLO γ rc with that of ∆γ rc ss , which effectively corresponds to using the same matching function Eq. (3.27) but without the last term.

Resummed DGLAP evolution matched to NNLO
As already discussed in Sect. 3, the ABF construction of the resummation of the anomalous dimension γ + relies on a double-leading (DL) part, which is based on the duality between the DGLAP and BFKL kernels (at the core of the resummation), and on a running-coupling (RC) part, which includes a class of subleading but very important effects which change the nature of the small-x singularity. The DL resummation is naturally performed at LO+LL or at NLO+NLL, which are obtained by combining together the (N)LO DGLAP anomalous dimension and the (N)LO BFKL kernel. Therefore, previous results on small-x resummation have always been presented at these orders.
As already mentioned in the introduction, it is of great importance being able to match the resummed result to fixed NNLO in order to obtain state-of-the art theoretical predictions. Matching the resummation to NNLO is in principle straightforward: starting from the NLO+NLL resummed result, one just needs to subtract its α s expansion up to O(α 3 s ), and replace it with the exact NNLO expression. While subtracting the NLO from the NLO+NLL is trivial, further subtracting the O(α 3 s ) term is not, due to the fact that the DL resummation is expressed in terms of implicit equations, which are usually solved numerically. One could think of different alternatives. One possibility is to expand the resummed result numerically, which however does not seem to be a reasonable option, as the numerical solution of the implicit equations is already challenging and slow, and one cannot hope in general to obtain sufficient precision in a reasonable amount of time from numerical techniques (unless further numerical developments are made 10 ). A second option is to construct a DL result starting from NNLO DGLAP and NLO BFKL, so that the result would be naturally NNLO+NLL. This option is itself non-trivial, as it requires the computation of a new class of double-counting terms between the two kernels, and has the undesirable disadvantage that the resummed result one obtains would differ by subleading NNLL terms from the original NLO+NLL.
In this work we have opted for a third, and perhaps more natural, option, namely expanding the resummed result analytically. Despite the rather technical nature of this computation, we find it illustrative to give its details in the following Sect. 4.1. Indeed, for instance, this exercise allowed us to find a small mistake in the original ABF construction of the DL part [37], which we also have inherited in our previous work [61], and which we correct here. Then, in Sect. 4.2, we present all the final expressions for the resummed splitting functions, providing a detailed explanation of the implementation of small-x resummation that constitutes the backbone of HELL, version 2.0.

Expansion of the Double Leading anomalous dimension
In the ABF construction, the DL resummed anomalous dimension γ DL , Eq. (3.1), is obtained from an implicit equation of the form where the function χ Σ (M, N, α s ) is a so-called off-shell BFKL kernel [37,61]. The DL anomalous dimension obtained through Eq. (4.1) assumes fixed coupling α s , and it thus receives a correction, Eq. (3.19), due to running coupling effects. This correction starts at O(α 4 s ) and it is therefore of no interest for the expansion of the DL result to O(α 3 s ). In the following, we explain how to construct a perturbative expansion of γ DL (N, α s ) defined by the implicit equation (4.1), focussing first on the simpler LL case, and moving next to the NLL case.

Expansion of the LL resummed result
We start from LL resummation for simplicity. We seek its expansion to O(α 2 s ), which would be needed to match LL resummation to NLO. The off-shell kernel at LO, needed for LL resummation, is given by [37] where the function χ s (α s /M ) is defined as the dual to the LO anomalous dimension γ 0 , the functionχ is the off-shell extension of the LO BFKL kernel with double counting with χ s subtracted, and the function Note that the LO anomalous dimension γ 0 (N ) that we use for the definition of χ s does not necessarily need to be the exact LO anomalous dimension. In fact, it can be replaced with an approximate expression with the same qualitative features and which preserves its small-x behaviour. This was already done in both Refs. [38,61], in slightly different ways, to cure a problem due to a branch-cut present in the n f = 0 case. Here, we adopt another, simpler, approximation, which circumvents the same problem and also solves another issue. Additionally, it allows us to exploit duality relations analytically, which is a great advantage from the numerical implementation point of view. Further details are given in App. B.3. In order to obtain the coefficient of the α s -expansion of the DL anomalous dimension, we substitute the formal expansion where χ 01 = χ s (0) = C A /π, from which it immediately follows Note that the O(α 0 s ) term cancels automatically, which confirms that indeed the LO part of γ DL-LO is given by the same γ 0 appearing in Eq. (4.3). Now, it happens that, due to the explicit form of and hence we findγ This might come as a surprise, however it does not. Indeed, the LL pole of the exact NLO γ 1 is accidentally zero, so the only part which is supposed to be predicted correctly by this kernel had to be zero. In principle there could be non-zero subleading corrections, which in practice are absent (at DL level -RC contributions do produce extra terms, see Eq. (3.17)). If we wish to match LL resummation to NNLO, we should expand to one extra order, but we are not interested in doing so, thus we move to the next logarithmic order.

Expansion of the NLL resummed result
For NLL resummation we need the NLO off-shell kernel   [37]), which did not contain the second term. In fact, the second term was not really necessary, as it is subleading, but then the argument of ψ 1 in the first term should be 1 − M + N , as the 2 comes from the subtraction of double-counting with the second term. In practice, however, we have verified that neglecting the second term and correcting the argument of the first leads to a kernel which is unstable close to the anticollinear pole M = 1, instability which is cured (resummed) by including the second term. We verified that the overall effect of this correction is mild, but not negligible. Finally, χ NLO mom (N, α s ) restores momentum conservation, in the same form as Eq.
The coefficients of the expansion of χ s,NLO are given by (see App. B.3) Now, from the definition ofχ 1 (see Ref. [61]) we immediately findχ which is N -independent, and thus also equal to the momentum conservation subtractionχ 1 (0, 1). Its value is (from Eq. (A.29) of Ref. [61]) On the other hand we have which is again N -independent. We can then rewrite Eq. (4.20) as This represents the final result for our expansion of the DL anomalous dimension. As a cross-check, we can now expandγ 2 about N = 0. Given that we have that up to NLL the singular behaviour ofγ 2 in N = 0 is given by (using f mom (0) = 0) which indeed reproduces the correct NLL pole of the known three-loop anomalous dimension γ 2 (N ) [101], while the LL 1/N 3 pole is accidentally zero. We stress that without the correction of the error in Eq. (4.16) the constant in Eq. (4.25) would change, and thus the NLL singularity of γ 2 (N ) would not be reproduced.

Resummed splitting functions matched to NNLO
In the previous sections we have computed the expansion of the DL result, and in Sect. 3 we have presented a new running coupling resummation and provided its α s expansion. We are now ready to construct the final expressions for the resummed anomalous dimension γ + , and write their expansions. With those, we can then also construct the resummed expressions of all the entries of the singlet anomalous dimension matrix in the physical basis [61], which by Mellin inversion give the singlet splitting functions.

Anomalous dimensions
As a first step, we need to add the running coupling contribution to the DL result, with the proper matching functions to cure the singularity mismatches. Note that adding the RC resummed functions ∆ DL-(N)LO γ rc (N, α s ) to the DL results violates momentum, which is further violated in the LO+LL by the matching function and in the NLO+NLL by ∆γ rc ss and its matching function. Momentum conservation can be restored by simply adding a function proportional to f mom (N ), Eq. (4.5). In summary, we have 12 where the various functions have been introduced in Sect. 3 and Sect. 4.1. Using the results in there, these expressions admit the following α s expansions and similarly for the LL result. In this way, the resummed anomalous dimension matched to the (exact) fixed order is given by where γ N n LO + (N, α s ) is the exact N n LO anomalous dimension. In Ref. [61] we only considered the "natural" contributions With the results of Eqs. (4.31) we can now also compute which are the resummed contributions needed for NLO+LL and NNLO+NLL. The previous equations are the primary ingredients which allow us to match NLL resummation to NNLO. Having the expansion of the eigenvalue anomalous dimension γ + , we can now construct the expansions for all the entries of the anomalous dimension matrix in the physical basis. For achieving this, the secondary ingredient is the expansion of the resummed γ qg entry of the evolution matrix. We refer the Reader to Ref. [61] for all the details on its resummation, and we report here its expansion in powers of α s , Note that in Ref. [61] two forms for γ NLL qg , differing by subleading terms, were considered and used to estimate an uncertainty of the resummation. Up to the order of Eq. (4.36), both expressions give identical results, the difference starting at O(α 4 s ). Using Eq. (4.36) it is possible to construct the resummed contributions (4.39) which are needed for NLO+NLL and NNLO+NLL resummations, respectively. The resummed contributions for other entries of the anomalous dimension matrix can be constructed in terms of ∆ n γ NLL + and ∆ n γ NLL qg , as described in Ref. [61].

Splitting functions
From the resummed anomalous dimension we can obtain resummed contributions for the splitting function matrix by Mellin inversion. Additionally, in order to ensure a smooth matching onto the fixed-order at large x, a damping is applied. Furthermore, we enforce exact momentum conservation on our final results by requiring the first moments of P gg + P qg and P gq + P qq to vanish. The final expressions are given by  4 . Note that, with respect to Ref. [61], we have introduced a further damping function, (1 − √ x) 4 , to ensure a smoother matching with the fixed-order at large x. From a numerical point of view, we proceed as follows: • With a new private code (available upon request) we produce the resummed anomalous dimension γ + . Specifically, we output ∆ 2 γ LL ( ) + and ∆ 3 γ NLL + , i.e. the ones with the expansion terms computed in this work subtracted, along the inverse Mellin integration contour, for a grid of values of α s and n f . This ensures that these contributions start at O(α 3 s ) and O(α 4 s ) respectively, as the subtraction is performed within the same code (subtracting later in a different code is of course possible, but subject to a higher chance of introducing bugs or numerical instabilities).
• The output of the first code (in the form of publicly available tables) is then read by the public code HELL, version 2.0 onwards. This code essentially computes the resummation of coefficient functions (and equally of the qg anomalous dimension) as described in Ref. [61], and partly summarized in Sect. 2.2. The splitting function matrix is then constructed and momentum conservation imposed. The objects which are computed are ∆ 2 P LL ig , ∆ 3 P NLL ig (i = g, q), ∆ 2 K hg , ∆ 2 c a,g , ∆ 2 c a,g (m) (a = 2, L) and ∆ 2 c CC a,g (m) (a = 2, L, 3) 13 on a n f , α s , x grid. Coefficient functions for additional processes will be also added in the future in the same form.
• These grids (again publicly available) are then read by the public code HELL-x, version 2.0 onwards, where the α s , x grid is interpolated (cubicly), the quark components of the splitting, matching and coefficient functions are computed by colour-charge relations, and ∆ 1 P LL ij , ∆ 2 P NLL ij , ∆ 1 K hg , ∆ 1 c a,i and ∆ 1 c (CC) a,i (m) are also constructed by adding the respective lower orders directly in x-space. For this, we need the analytic inverse Mellin transforms of the non-trivial expansion terms in Eq. (4.31a), (4.31b) and (4.31c), as well as the analogous for the coefficient functions. The latter was already done in the massless case in Ref. [61]; explicit results for the massive case are presented in App. A.3. Explicit results for the splitting functions are presented in App. B.4.
We underline that the second step is the slowest, as the HELL code needs to compute several integrals for the varius grid points and the various functions to be resummed. The last step, performed by the HELL-x code, is instead extremely fast, as it simply amounts to an interpolation and the evaluation of simple functions. For practical applications of our results, it is sufficient to use the HELL-x code, making the inclusion of small-x resummation very handy. As far as PDF evolution and construction of DIS observables is concerned, the code HELL-x has been included in the APFEL code [91], which can be used to access all our results and construct resummed predictions for physical observables.

Numerical results
In order to illustrate the capabilities of HELL 2.0, we present here some representative results for the small-x resummation of splitting functions and DIS coefficient functions obtained with the 13 Massive DIS coefficient functions are further sampled for various values of the quark massess. The results also include an uncertainty band, as described in the text. The plots are for αs = 0.2 and n f = 4 in the Q0MS scheme. We note that difference between Q0MS and MS for the fixed-order results is immaterial at this accuracy.
techniques described in Ref. [61] and improved as described in the previous sections. Moreover, we also show new results for the coefficient functions with massive quarks.

Splitting functions
Let us start with DGLAP evolution. With respect to our previous work [61] we have made substantial changes in the resummation of the anomalous dimensions, mostly due to the treatment of running coupling effects, as described in Sect. 3. Additionally, we are now able to match the NLL resummation of the splitting functions to their fixed-order expressions up to NNLO, as presented in Sect. 4. In Fig. 2 we show the fixed-order splitting functions at LO (black dotted), NLO (black dashed) and NNLO (black dot-dot-dashed) compared to resummed results at LO+LL (green dotted), NLO+NLL (purple dashed) and NNLO+NLL (blue dot-dot-dashed). In principle, we also have the technology for matching LL resummation to NLO, but this is of very limited interest, so we do not show these results here (they can be obtained from the HELL-x code). The gluon splitting functions P gg and P gq are shown in the upper plots, and the quark ones P qg and P qq are shown in the lower plots (the latter two start at NLL so the LO+LL curve is absent there). All splitting functions are multiplied by x for a clearer visualization. The scheme of the resummed splitting functions is Q 0 MS (the fixed-order ones are the same in both MS and Q 0 MS at these orders). The number of active flavours is n f = 4, and the value of the strong coupling is α s = 0.2, corresponding to Q ∼ 6 GeV. Note that for such value of Q in a VFNS one usually has n f = 5 active flavours; however, the difference between the results in the n f = 4 and n f = 5 schemes at the same value of α s is modest, and our choice allows to directly compare with previous results presented in the literature.
The results of Fig. 2 can be compared directly to the ones presented in our previous paper [61]. It can be noticed that the LO+LL result is rather different: in our new implementation this curve is lower than in the previous version. This is entirely due to the new treatment of running coupling effects, which clearly differs by subleading logarithmic terms. At the next order, NLO+NLL, there is again a difference with respect to our previous work. As before, this is due to subleading terms, which are now NNLL, and so, as expected, they lead to smaller discrepancies. Indeed, NLO+NLL results are much closer to the ones of our previous implementation. We recall that the new version of the NLO+NLL results also includes the correction of an error in the original expressions [37], as detailed in Sect. 4.1.2, which has a non-negligible impact on the result, even though the effect is not as large as the one induced by the new treatment of running coupling effects.
The notable novelty is the presence of the NNLO+NLL curve. The asymptotic small-x behaviour is identical to the NLO+NLL curve, except for a constant shift, which represents a term of the form α 3 s /x in the splitting functions. Indeed, this term is NNLL, and it was therefore not correctly captured by the NLO+NLL result. Its impact is larger in the gluon splitting functions P gg and P gq , while it is rather small for the quark splitting functions P qg and P qq . At larger x, the NNLO+NLL curves smoothly match onto the NNLO result. For P gg and P gq this happens already at x ∼ 10 −2 . This is due to the fact that the dip at x ∼ 10 −3 ÷10 −4 , which is a known feature of the resummed result at moderate x (see e.g. [30,31,37,102]), is determined by the NNLO logarithmic term, which goes down and dominates at moderate values of x, before the onset of the smaller x asymptotic behaviour, which goes up. Hence, when matching to NNLO, the initial descent of the splitting functions is automatically described, and the resummed result deviates only when the rise due to the asymptotic small-x behaviour sets in. Note that, because of this difference between NLO+NLL and NNLO+NLL, we expect the latter to be much more accurate, especially in the region x 10 −5 , where the majority of DIS data lie.
The resummed curves are supplemented with an uncertainty, which aims to estimate the size of subleading logarithmic effects. As far as P qg and P qq are concerned, it is defined exactly as in Ref. [61], namely symmetrizing the difference between our default construction of the γ qg resummation which uses the evolution function in Eq. (B.1), and what we obtain by switching to a simpler and formally equally valid version, Eq. (B.3). The uncertainty band is bigger than in our previous work just because the LL anomalous dimension used to resum γ qg differs in the treatment of running coupling effects. Because of the way P gg and P gq are constructed, Eq. (4.41), these splitting functions inherit the uncertainty band of P qg at NLL. However, because the bulk of resummed contribution to these entries comes from γ + , we have decided to also account for the uncertainty due to subleading contributions to γ + . This was not considered in our previous work. In order to construct this band, we use the same kind of variation used for γ qg . Specifically, we symmetrize the difference between the result obtained using Eq. (3.3) for the resummation of running coupling effects (our default) and the variant obtained using the simpler, yet equally valid Eq. (3.4). The uncertainty bands from both sources are then combined in quadrature. At LL, there is no contribution from γ qg , and the whole resummed curve is given by γ + : in this case the uncertainty band is just determined from the variation in the construction of the latter.
We note that there is nice overlapping between NLO+NLL and NNLO+NLL bands for the P qg and P qq splitting functions, giving us a good confidence that they appropriately represent the uncertainty from missing subleading logarithmic orders. In contrast, the uncertainty band on P gg and P gq does not fully cover the effects of higher orders in the initial small-x region, 10 −4 x 10 −2 , as demonstrated by the fact that NLO+NLL and NNLO+NLL do not overlap there. However, this effect is mostly driven by the largish NNLL effects at O(α 3 s ), which are those that are included in the NNLO+NLL but not in the NLO+NLL results. At higher orders the effects of subleading logs in this region is likely to be smaller. In support of this hypothesis, we can note that the distance between NNLO+NLL and NNLO for x ∼ 10 −2 is significantly smaller than the distance between NLO+NLL and NLO, in the same region. Thus, we believe that, while the uncertainty on the NLO+NLL result is not satisfactory in the intermediate x region, the uncertainty on the NNLO+NLL should be reliable.
These plots are also instructive to study the stability of the perturbative expansion. By looking at the fixed-order splitting functions, we see that small-x logarithms start being dominant already at x 10 −2 , where the logarithmic term of the NNLO contribution sets in. We note that the small-x growth could have been in principle much stronger. Indeed, the leading logarithmic contributions have vanishing coefficients both at NLO and NNLO and the sharp rise of the NNLO splitting function is driven by its NLL contribution α 3 s x −1 log x. These accidental zeros are not present beyond NNLO and so we expect the yet-unknown N 3 LO splitting functions to significantly deteriorate the stability of the perturbative expansion because of their α 4 s x −1 log 3 x growth at small x. Therefore, we anticipate that the inclusion of the resummation to stabilize the small-x region will be even more crucial at N 3 LO.

DIS coefficient functions
We now move to DIS coefficient functions and we first present updated predictions for the massless coefficients, i.e. assuming that there are only n f massless quarks and no heavy quarks. The construction did not change compared to our previous work [61], but the input LL anomalous dimension used for computing the resummed coefficients did, as explained in Sect. 4. 14 The updated results are shown in Fig. 3 for C L,g (left) and C 2,g (right). The quark contributions at small x are very similar (due to colour-charge relation) and are not shown. We observe some differences with  Figure 4. Same as Fig. 3, but for the massive coefficient functions cL,g (left) and c2,g (right), for both charm production (upper plots) and bottom production (lower plots) in NC. The charm production plots are for n f = 3 and αs = 0.28, corresponding to Q ∼ 2 GeV, slightly above the charm mass mc = 1.5 GeV.
The bottom production plots are for n f = 4 and αs = 0.2, corresponding to Q ∼ 6 GeV, slightly above the bottom mass m b = 4.5 GeV.
respect our previous work, although within uncertainties, for the coefficient function C 2,g due to the modified running-coupling resummation. These changes appear to have instead a small numerical effect on C L,g . The other noticeable difference with respect to our previous results is the size of the theoretical uncertainty, which is now larger: this effect is entirely due to the different LL used in the construction, and is therefore ultimately due to the treatment of running-coupling effects. We now move to the new results which include mass dependence. We first show in Fig. 4 the analogous of Fig. 3 for the massive unsubtracted coefficient functions, both for charm production and for bottom production close to the production threshold. As usual in theory papers, we define these contributions as the ones for which the heavy quark is struck by the photon (at these energies, the Z contribution in NC and the CC production mechanism are negligible). We call generically these contributions c a,i , with a = 2, L and i = g, q, of which the functions ∆ c a,i defined in Sect. 2.1.1 are the resummed contributions. For charm production (upper plots) we use α s = 0.28, corresponding to Q ∼ 2 GeV, which is a scale right above the charm mass assumed to be m c = 1.5 GeV, while for bottom production we use α s = 0.20, corresponding to Q ∼ 6 GeV, right above the bottom mass assumed to be m b = 4.5 GeV. The number of active flavours is set to be n f = 3 for charm production and n f = 4 for bottom production, i.e. the massive quark is treated as heavy and its collinear logarithms are not factorized. In particular, the massive coefficients for bottom production are those contributions which should be added to the corresponding massless coefficients in the same charm production -CC -x = 10 -3 Figure 5. Resummed contributions to the coefficient functions including mass effects in neutral current, ∆2cL,g(mc) and ∆2c2,g(mc) (upper plots), and in charged current, ∆2c CC L,g (mc), ∆2c CC 2,g (mc) and ∆2c CC 3,g (mc) (lower plots), shown as a function of Q from mc to 100 GeV. The left plots show the results for x = 10 −6 , the right plots for x = 10 −3 . The value of n f changes from 4 to 5 at the bottom matching scale taken to be equal to m b . The massless result, to which the massive coefficients tend at large Q, is also shown as dotted curves (except for ∆2c CC 3,g (mc), which tends to zero). n f = 4 scheme, Fig. 3, which instead assumed only coupling to light quarks, to obtain a complete prediction (see e.g. the decomposition Eq. (2.16) at resummed level). We observe that the effect of adding the bottom production contribution to the purely massless contributions is a rather small correction for F L , while it is comparable in size to each individual massless contribution for F 2 . The pattern observed in Fig. 4 between fixed-order and resummed contributions is very similar to that of the massless results in Fig. 3. The most notable difference is the visibly larger effect of resummation for charm production, accompanied by a larger uncertainty band. This effect is entirely due to the smaller scale, i.e. the larger value of α s , used in the charm production plots. Another interesting feature of these massive coefficient functions is the very visible presence of the physical threshold for heavy quark production, which lies at x = x th ≡ 1/(1 + 4m 2 /Q 2 ). Because of our choice of scales, x th ∼ 0.3 for both processes.
The results presented so far do not include the resummation of collinear logarithms due to massive quarks. For the scales considered, which are in both cases larger than the heavy quark mass, these collinear logarithms are already usually resummed in most implementations of VFNSs. We now thus consider the scenario in which a VFNS is used and heavy-quark collinear logarithms are resummed. Since at fixed order there are various incarnations of VFNSs, differing just by subleading effects but nonetheless being practically different (see e.g. discussion in Sect. 2.2.3), we prefer to focus on the resummed contributions only. We focus on the charm-production case, with m c = 1.5 GeV. For the sake of this study, we find more interesting to show a plot as a function of the momentum transfer Q, in order to emphasize the importance of mass effects at different scales. Therefore, in Fig. 5 we plot the resummed contributions ∆ 2 c 2,g (m c ) and ∆ 2 c L,g (m c ) in NC and ∆ 2 c CC L,g (m c ), ∆ 2 c CC 2,g (m c ) and ∆ 2 c CC 3,g (m c ) in CC 15 as a function of Q for two representative values of x, namely x = 10 −6 (small) and x = 10 −3 (moderate). The plot starts from Q = m c , where n f = 4, and then it transitions to n f = 5 when crossing the bottom threshold (assumed to be at m b = 4.5 GeV). At the transition point a small discontinuity appears, due to the different value of n f used in the computation of the LL anomalous dimension. This discontinuity is a standard consequence of the scheme change, and does not constitute any practical problem in the computation of physical observables.
At large Q, the massive resummed coefficient functions (which are the collinear subtracted ones) tend to the massless results, shown in dotted style. It is clearly visible that charm mass effects are significant for small Q (10 ÷ 30) GeV, and are more pronounced at larger x, where however the effect of resummation is smaller. Charm mass effects are also stronger in the NC case than in the CC case. In practice, massive corrections on the resummed coefficient functions are a small effect on the full structure function, especially when resummation is matched to NNLO. Still, keeping into account these mass effects is important for an accurate description of the low-Q data, and in particular for the charm structure function F c a , a = 2, L, which is entirely determined by the charm coefficient function.
In the upper plots of Fig. 5 we are showing the full NC coefficients, namely the sum of the contributions from photon, Z and photon-Z interference to the structure functions, normalized to the photon couplings. It is interesting to investigate how much the various terms contribute to the full result. To do so, we show in Fig. 6 the ratio of the individual contributions to the resummed contribution to the structure function F 2 , ∆ 2 c 2,g (m c ). We stress that if the axial contribution proportional to g 2 A which remains after factoring out the g 2 V + g 2 A coupling to the Z were absent, then the resummed coefficient function for the various contributions would be identical, up to the overall coupling, and the ratio of the various contributions would be independent of x and of the observable. However, since this axial contribution is non-zero, a small dependence remain. However, for x 10 −3 , the differences between F 2 and F L are very small, and reducing with lowering x. Thus, the plot in Fig. 6, obtained from F 2 for x = 10 −4 , remains pretty much unchanged for lower x and for F L . From the plot we see that the contribution from the Z boson is basically negligible for all Q 10 GeV, and becomes of the same size of the photon contribution at the Z peak. The photon-Z interference dominates over the pure Z-exchange contribution below the Z peak. The axial contribution proportional to g 2 A , which is a new result of our computation, turns out to be mostly insignificant, as it gives a small contribution (a few percent at most) for scales Q where small-x resummation is further suppressed by a smaller strong coupling.

Conclusions
In this paper we have performed a comprehensive study of the high-energy, i.e. small-x, resummation in deep-inelastic scattering of a lepton off a proton. In particular, we have collected all the ingredients to perform NLL resummation matched to fixed-order up to NNLO.
In order to achieve this we have considered the resummation of splitting functions, which govern DGLAP evolution of parton densities. With respect to our previous work, we have modified the way running coupling corrections are treated and we have managed to match the resummation to NNLO, thus obtaining state-of-the art NNLO+NLL results for the splitting kernels. While fixedorder predictions at NNLO exhibit instabilities at small x due to large logarithms, the resummed results are stable and at small x appear to be much closer to NLO than NNLO. Furthermore, our results can be easily extended to match NLL to fixed N 3 LO, when it becomes available. In this case we expect the resummation to have an even more substantial effect because of the larger fixed-order instabilities at small x appearing at this order.
We have also considered the resummation of DIS partonic coefficient functions. In order to obtain reliable results in a wide range of x and Q 2 we have studied small-x resummation in the context of a variable flavour number scheme in which heavy and light quark coefficient functions are matched together. In this context, we have considered mass effects originating from both the charm and the bottom quarks. We have produced NNLO+NLL results for the coefficient functions relevant for F 2 and F L for neutral-current DIS, considering the effect of both a virtual photon or a Z boson exchange, as well as charged-current processes. If all quarks are massless, the structure function F 3 is purely non-singlet and therefore is not enhanced at small-x. However, we have found that in the charged-current case with W boson exchange, if at least one of the quarks interacting with the W is massive there is a non-zero contribution at small x to the parity violating structure function F 3 . We have also noted that in neutral-current DIS with massive quarks there is a difference between the γ exchange or Zγ interference and the pure Z exchange, other than the overall coupling.
We have implemented all these new results in a new version of our code, HELL version 2.0. A fast interface to these results is available through a new version of its companion code, HELL-x version 2.0. Both codes are publicly available at the webpage www.ge.infn.it/∼bonvini/hell. The fast HELL-x 2.0 code can be directly used to compute PDF evolution and DIS cross sections through the public code APFEL [91].
The main motivation behind this work was to compute and implement all the ingredients that are necessary to perform a state-of-the-art fit of parton distribution functions which consistently includes small-x resummation both in the evolution of the parton densities and in the coefficient functions. This task is being now pursued by the NNPDF collaboration. Preliminary and very encouraging results have been presented in [66] in the case of a PDF fit that includes DIS-only data. For the near future, we look forward to implementing other processes in HELL, the most relevant of which in the context of PDF extractions is the production of lepton pair via the Drell-Yan mechanism. We conclude by noting that the resummed results produced by HELL, both splitting functions and coefficient functions, are supplemented by a band representing the theoretical uncertainty due to missing higher-logarithmic orders. This information can (and should) be used in phenomenological studies and it could be also be fed into PDF fits together with other sources of theoretical uncertainty. However, the debate about how to best include theoretical uncertainties into PDF fits is not settled yet, see for instance [103][104][105][106]. Figure 7. Diagrams entering the computation of the DIS coefficient function with off-shell incoming gluons at lowest order. The (off-shell) vector boson can either be a photon, a Z or a W . In the latter case, the quark changes flavour after interacting with it, and thus does the mass.
where g s , g V , g A are the strong, vector and axial couplings respectively, and t = (k − p 3 ) For photon-induced DIS, the vector coupling is just the quark electric charge g γ V = e Q and the axial coupling is zero. For Z-induced DIS, the vector and axial couplings depend on the quark isospin and are given by where θ w is the weak mixing angle. For W -induced DIS the vector and axial couplings depend on both quark flavours and are given by being V ij the CKM matrix. Note that we are assuming that the vector boson is just V , and so the computation will not cover explicitly the photon-Z interference: however, this case is easily obtained in the final results by simply replacing g 2 V → 2g Z V g γ V and g 2 A → 0 (the combination g V g A appears only in F 3 , the gluonic contribution of which is zero in neutral current).
In the high-energy regime we are interested in, we decide to parametrize the kinematics of the process in terms of dimensionless variables z 1 , z 2 , τ ,τ and transverse vectors 16 q ⊥ , k ⊥ and ∆ ⊥ defined by where p 1 = √ s 2 (1, 0, 0, 1) and p 2 = √ s 2 (1, 0, 0, −1) are light-cone vectors. We also define the momentum transferred Q by q 2 = −Q 2 . It is important to note that this definition is valid only in the high-energy limit, where s M 2 V .
We define the parton-level off-shell hadronic tensor as the squared of the matrix element Eq. (A.2) averaged over the off-shell gluon polarizations [45,51,107] This object contains all the information on the DIS cross-section from the hadronic side. It can be decomposed into different contributions with a given tensor structure as 17 The various contributions to the structure functions F 1 , F 2 and F 3 (and F L = F 2 − 2xF 1 ) are obtained from the respective counterparts in the hadronic tensor W 1 , W 2 and W 3 . These, in turn, can be obtained from the full tensor W µν using suitable projector operators. To this end, we define where |A|  17 Note that in the photon-mediated DIS the contributions W 4 and W 5 can be related to W 1 and W 2 through Ward identities. However, since we want to cover also the more general Z-and W -mediated DIS processes, we must keep them separate. Nevertheless, their contribution to the DIS cross section, as well as the one from W 6 , is of the order of the lepton mass and thus negligible, and will not be further considered.
These squared amplitudes must be integrated over the final state two-body phase space. In terms of the kinematic variables Eqs. (A.5) we can express it as where we have further defined We recall that bold symbols represent the two-dimensional components of a transverse vector. The partonic off-shell coefficient functions (in Mellin space) for the three structure functions we are interested in are given in terms of the squared amplitudes of Eqs. (A.8) by where we have introduced the variable η = Q 2 ν (A.14) and we are expressing the result in terms of dimensionless ratios Note that, as explained in the text, we are only interested in the N = 0 Mellin moment.
To carry out the integrations in Eqs. (A.13) it is useful to express the following combinations in terms of the phase-space variables. In addition, it is convenient to use the following Feynman parametrizations In this way, most integrations are easy to perform, 18 leaving the results in the form of a double integral in y and τ . In order to simplify this integration (and to make contact with previous literature) it is also convenient to perform the change of variables and express the results as integrals over x 1 and x 2 . General results for C a (0, ξ, ξ m1 , ξ m2 ), a = 2, L, 3 are rather long and will not be reported here; in the next section we focus on the physical combinations that are relevant for neutral and charged currents, where the masses are either equal or at least one of them is vanishing.

A.2 Results
In this section we collect the results of the off-shell coefficient functions for neutral and charged currents, contributing to the three structure functions F 2 , F L and F 3 . The fully massless case m 1 = m 2 = 0, which is common to neutral and charged currents, is of course the simplest limit and yields These results coincide with those presented in Ref. [49], even though the longitudinal coefficient function was not written explicitly there. An equivalent (but simpler) integral form for C L (0, ξ, 0, 0) was also given in Ref. [61]. We now move to the case in which both masses are equal, m 1 = m 2 ≡ m, which is relevant for neutral current. The results read Note the presence in the above expressions of a term proportional to g 2 A /(g 2 V +g 2 A ). This contribution is present only when the vector boson is a Z, while for photon and photon-Z interference this term is zero. This axial contribution is a new result. The remaining of the expressions were already known [44,45,90], however an explicit integral form of this kind for C L is presented here for the first time. Of course, the massless ξ m → 0 limit of Eqs. (A.19) reduces to Eqs. (A.18).
Finally, we move to the case in which one quark is massless (say, m 2 = 0) and the other is massive (m 1 ≡ m), which is relevant for charged current. According to the definition of the process, Eq. (A.1), this choice corresponds to the production of a heavy anti-quark. Here, it is most convenient to leave integration over τ untouched and to change variable only from y to x 1 . We thus obtain To the best of our knowledge, these results are all new. Most notably, C 3 does not vanish in this case. Note that choosing m 2 = m, m 1 = 0 corresponds to charge-conjugating the final state, thus producing a heavy quark. Therefore, C a (0, ξ, 0, ξ m ) = C a (0, ξ, ξ m , 0) for a = 2, L, while there is a sign change in the parity-violating coefficient, C 3 (0, ξ, 0, ξ m ) = −C 3 (0, ξ, ξ m , 0) (see also Ref. [108,109] We now consider the Mellin transform with respect to ξ of these results. This is particularly useful for studying the massless limit of the resummed result, and for asymptotic expansions. We denote the Mellin transform with a tilde, and replace the argument ξ with the Mellin moment M : (N, ξ, ξ m1 , ξ m2 ), a = 2, L, 3. (A. 21) In the massless case we find which reproduce the results of Ref. [49]. In the massive case in neutral current we obtaiñ Apart from the contribution proportional to the axial coupling g 2 A , which is new, the other terms reproduce the results of Refs. [44,45,90]. Finally, the massive case in charged current is given bỹ

A.3 Special limits
We find instructive to study the above expressions in two particular limits, namely the massless limit and the limit M → 0. The latter is useful to construct the fixed-order expansion of the resummed result.

A.3.1 Massless limit
We can use the Mellin forms to study the massless limit ξ m → 0 of the resummed expressions. We have for neutral current , can be obtained in the following way. We first split the function as the product of three different factors which we write in integral form: Then we change variable from t to ξ = 4ξ m ty/x to write the product in the form of a Mellin transfrom:K At this point we can simply read off the inverse Mellin transform in integral form 34) which has been computed explicitly in the second line. The ξ-derivative of this expression is Eq. (2.42).

A.3.2 Fixed-order expansion and collinear singularities
The Mellin form of the off-shell coefficients can be also used to compute expansions in M = 0, which are needed for computing the α s expansion of the resummed results. For matching up to NNLO, i.e. O(α 2 s ), we need the expansions up to O(M ). In the massless case we havẽ while in the massive case in neutral current we obtaiñ having defined the harmonic polylogarithm In the charged-current case, with only one massive quark, we have instead the following expansion: (A.38c) In each of the above equation, the pole in M = 0, where present, identifies the collinear singularity and is given by the LO P qg in Mellin space times the LO non-singlet process q + W → q . The latter has non-trivial mass dependence in charged current, Eqs. (A.38), since in this case the final state quark q is massive. Note that in this case even F L has a non-vanishing contribution at LO, which is proportional to the mass and thus vanishes in the massless limit. The O(M 0 ) terms in these expansions, after subtraction of massless collinear singularities as in Eqs. (2.31a) and (2.35), reproduce the known O(α s ) contributions [110,111]. Higher-order corrections exist both for neutral and charged currents (see e.g. [112,113]), however not in a form we could compare our expansion to.

B.1 The evolution function U
A key ingredient of the formalism for resumming coefficient functions is the evolution function U (N, ξ), defined in Eq. (2.32). As discussed in Ref. [61], computing it exactly with the resummed anomalous dimension γ + (specifically, with the LL anomalous dimension) requires integrating γ + over all values of α s from 0 to ∞, which is numerically inconvenient. Therefore, following Refs. [38,50], we use an approximate expression for the evolution factor, where the anomalous dimension is assumed to depend on α s only at LO, with 1-loop running. This leads to the ABF evolution factor [61] U ABF (N, ξ) = 1 + r(N, α s ) log ξ γ+(N,αs)/r(N,αs) Note that the ratio r(N, α s ) is such that the approximation reproduces the correct derivative of γ + in α s . However, this effect is strictly speaking beyond the formal accuracy we work with, so one could ignore it and replace This variant is used to construct the uncertainty band on our resummed predictions. As a simpler approximation we could also consider the fixed-coupling limit, in which all the scale dependence in α s is ignored and the evolution factor becomes simply In this case our formula for computing the resummation of coefficient functions simply reduces to a Mellin transformation with moment γ + . The integration range in the off-shellness ξ in the resummation formula extends to all accessible values between 0 and ∞. In the running-coupling case, α s is computed at ξQ 2 in U , Eq. (2.32), so at some small value of ξ the Landau pole is hit, and the integration must stop there. With 1-loop running (and also at higher loops if the expanded solution for the running coupling is used) the position of the Landau pole is given by and ξ-integration must be limited to the region ξ ≥ ξ 0 . In the fixed-coupling limit, α s is frozen at its value in Q 2 , so all values of ξ are in principle accessible, i.e. ξ 0 f.c.
= 0. In Ref. [61] we derived the resummation formula which had originally the form (schematically) Then, for convenience in the numerical implementation, we integrated by parts to get where the boundary term at infinity vanishes thanks to C, and the boundary term in ξ 0 is assumed to vanish because The latter assumption is not always true. However, at the leading logarithmic accuracy, which is the only accuracy on which we have control at the moment, the resummed result is only governed by the fixed-coupling anomalous dimension γ s , dual of the LO BFKL kernel. Thus, the formula Eq. (B.4) applies, with γ + = γ s . To obtain the resummed coefficient function in momentum space, an inverse Mellin transform has to be computed. This amounts to integrating in N over an imaginary contour with abscissa to the right of the small-x singularity, which in the case of γ s is placed in N = α s c 0 , with c 0 given in Eq. (3.18). Along such contour, the real part of γ s is always positive, and therefore U f.c. (N, ξ 0 = 0) = 0. Therefore, at the accuracy we are working with, the boundary term indeed vanishes.
In practice, however, we include in our resummation subleading contributions which spoil the condition Eq. (B.8). Indeed the anomalous dimension γ + that we use has a more complex structure than γ s . Additionally, we use the approximation Eq. (B.1), which is typically finite in ξ = ξ 0 :
Note that using the variant Eq. (B.3) U ABF (N, ξ 0 ) is either 0 or ∞ depending on the sign of the real part of γ + . This implies that the two formulations Eqs. (B.6) and (B.7) will give in general different results, due to the neglected non-zero boundary term. We have indeed verified this numerically. Despite the fact that this difference is subleading log, and hence either result is formally equally valid, this difference in the formulations is undesirable.
In this work we propose to solve the ambiguity by modifying the evolution function with suitable higher-twist terms such that we always have U (N, ξ 0 ) = 0. To do so, we use the evolution function It is easy to verify that the damping function D higher-twist (ξ) vanishes in ξ 0 and smoothly tends to 1 in ξ = 1, with all derivatives vanishing in ξ = 1. Moreover, it is clearly higher-twist, i.e. non analytical in the coupling α s , so it does not influence the perturbative expansion of the evolution factor (which is used for the matching of the resummed expressions to fixed-order). Using this new damped evolution function, we find that the results obtained using Eq. (B.7) are indistinguishable from those obtained with the undamped function, which confirms that the results of Ref. [61] are unaffected (from the point of view of U ). On the other hand, results obtained using Eq. (B.6) are now identical to those obtained using Eq. (B.7), as they must, since now the boundary term is identically zero by construction.
From the point of view of the numerical implementation, we observe that the N dependence of the resummed coefficient functions is all contained in U , Eqs. (B.6) or (B.7). Therefore, we can first compute the inverse Mellin transform of U (N, ξ) Eq. (B.10) as function of ξ, U (x, ξ), which we tabulate for various values of α s , x and ξ, and we subsequently use it to compute the ξ integration for each observable. In HELL v2.0 both the old methodology (which integrates first in ξ and then in N ) and the new one (integration order inverted) are implemented, and give of course identical results (within numerical integration errors). The new implementation is faster.

B.2 Implementation of kinematic theta functions at resummed level
In massive coefficient functions, kinematic constraints for the production of the massive final state are implemented through theta functions appearing in the coefficient functions. For instance, in the case of DIS neutral-current structure functions, the theta function has the form θ(X − x), with X = 1/(1+4m 2 /Q 2 ) and where x is the Mellin integration variable. 19 The very same theta function appears also in the off-shell coefficient, as it depends only on the kinematics of the final state. In the resummation procedure, the off-shell coefficient function is Mellin transformed with respect to x and then the result is evaluated in Mellin moment N = 0. The last step (which is strictly speaking not necessary) loses track of the theta function, and the inverse Mellin transform of the resummed on-shell coefficient is non-zero also in kinematically unaccessible regions.
A possible solution to restore the kinematic theta function is simply to avoid computing the off-shell coefficient in N = 0. This is possible, however, it is not convenient for at least two reasons: the first is that all expressions and calculations become significantly more complicated, and the second is that for consistency this should be done also in the massless case. The latter requirement is necessary in the construction of the collinearly resummed coefficient functions, otherwise the massless limit of the (collinear subtracted) resummed massive on-shell coefficients would not tend to the massless ones, Eq. (2.25).
Therefore, we seek a solution which restores the theta function in the resummed approach, while keeping using the off-shell coefficient in N = 0. The implementation must satisfy three requirements: • the theta function should be restored without affecting the logarithmic accuracy of the result; • the x → X limit must be smooth; • in the massless X → 1 limit the effect must disappear completely.
The first requirement is obvious. The second one is perhaps not mandatory, but it is satisfied in fixed-order results, and avoids sharp transitions between results. The latter requirement is instead needed for a correct implementation of the resummation of collinear logarithms at small-x resummed level.
We have investigated different options for the restoration of the theta function such that the requirements above are satisfied. We report here the two main alternatives that we consider, which act on N space and on x space, respectively. The N -space approach consists in multiplying the integrand of Eqs. As explicitly indicated, this term is manifestly subleading, and reproduces the theta function θ(X − x) thanks to the X N term. In the massless limit it reduces to Θ N (N, 1) = 1, as required. It can be also verified that, in full generality, the inverse Mellin transform of the resummed coefficient function obtained with this extra function vanishes smoothly as x → X. The alternative implementation in x space is obtained by multiplying the final resummed coefficient function in x space by the function The function in squared brackets ensures smooth x → X limit, and it is clearly subleading. In the massless limit X → 1 it reduces to Θ x (x, 1) = θ (1 − x), as required.
The two alternatives are formally equally valid, but lead in general to different numerical results. For practical reasons, we opt for the x-space implementation, Eq. (B.13). In this way restoring the theta function can be done at the very end, giving full flexibility for the implementation of the resummation. For instance, it is possible to precompute the inverse Mellin transform of the evolution function, U (x, ξ), as described in Sect. B.1, speeding up the computation of resummed massive coefficient functions. This would not be possible using the N -space formulation, Eq. (B.12), as in this case the N dependence of ξ-integrand of resummed coefficient functions would include the Θ N (N, X) term, so the Mellin inversion would not act on U (N, ξ) only.

B.3 A convenient approximate form for the fixed-order anomalous dimension
In the construction of the off-shell kernel for LL and NLL resummation, we need the dual of the fixed LO or NLO anomalous dimension, denoted χ s and χ s,NLO , respectively. These two functions provide the resummation of collinear (and anticollinear) contributions in the DL BFKL kernel. If the duals are computed from approximate expressions of the fixed-order anomalous dimensions, the resummation in the BFKL kernel is only approximate, and one cannot claim to have exact leading or next-to-leading logarithmic accuracy in BFKL. However, our goal is to reach LL and NLL accuracy in the resummed anomalous dimension, not in the BFKL kernel. For this, the BFKL kernel has just to be correct at fixed LO or NLO, since by duality this fully determines the LL and NLL contributions of the DGLAP anomalous dimension. The reason why also the BFKL kernel is being resummed is just that the resummation stabilizes its perturbative expansion, which is otherwise highly unstable close to the collinear and anticollinear poles. Therefore, from the point of view of the accuracy of the result, it is well possible to use an approximate expression for the LO and NLO anomalous dimensions, provided their LL and NLL parts are correct, as they correspond by duality to contributions of O(α s ) and O(α 2 s ) in the BFKL kernel, which need to be correct. Then, once the fixed-order part of the resummed (N)LO+(N)LL anomalous dimension is subtracted, the resummed contributions Eqs. (4.34a), (4.34b) can be added to the exact fixed-order anomalous dimension, restoring the correct fixed-order part.
This approach was already considered in both Refs. [38,61], with two different implementations. The basic motivation was that the exact fixed-order anomalous dimension γ + , being the eigenvalue of a matrix, contains a square-root branch-cut, which is inherited by the DL anomalous dimension and would give rise to a spurious oscillating behaviour when performing an inverse Mellin transformation. The approximate γ + implemented in Ref. [61], which is a somewhat simplified version of the one originally proposed in Ref. [38], was simply obtained by taking γ gg computed for n f = 0 (which is then also the eigenvalue, as there are no quarks and the matrix reduces to a single entry) and adding the n f -dependent contributions of the exact γ + restricted to LL and NLL. This procedure ensures that the resulting anomalous dimension reproduces the LL and NLL behaviour of the exact one, but it behaves as γ gg elsewhere in the N plane. This implies in particular that the anomalous dimension grows (negatively) as log N at large N . Note also that this construction violates momentum.
We observe that the large-N logarithmic growth of γ + is problematic. Indeed, the dual function χ s (or χ s,NLO ) grows exponentially for negative M as |M | gets larger. The DL kernel, by duality, should then be able to reproduce the logarithmic growth at large N . However, the DL kernel does not only contain χ s , but also the fixed-order BFKL kernel, which contains poles for all integer values of M . Therefore, by duality, the large N behaviour of the DL anomalous dimension is determined by the rightmost M pole for negative M , which is in M = −1, which implies that the DL anomalous dimension tends to −1 as N → ∞. This problem was ignored in previous works, as the M = −1 pole represents a higher twist contribution, and the practical effect is almost negligible. However, it would be ideal to avoid this issue. One option would be to act on the DL BFKL kernel, hacking it such that the poles for negative M are no longer present. While we have tried this solution, we think that it is not the best approach. A significantly better solution is obtained if the anomalous dimension used in the duality relation does not grow at large N , but rather it goes as a constant larger than −1, such that χ s never hits the rightmost negative M pole: in this way, the large-N behaviour of the DL anomalous dimension is determined by χ s itself, and hence corresponds to the one of the input anomalous dimension. Thus, here we propose a new approximation for the fixed-order anomalous dimension. We require that the LL and NLL behaviour is reproduced, that momentum conservation is preserved exactly, and that at large N it behaves as a constant greater than −1. Given that the LO and NLO anomalous dimensions behave close to N = 0 as Note that a 20 is formally NNLL, so it could in principle be ignored (similarly, for LL resummation one could ignore a 10 ); in practice, including both a 11  (here we generically use the name χ s for representing both the dual to the LO anomalous dimension and the dual to the NLO one, previously called χ s,NLO ). This represents a great advantage from the point of view of the numerical implementation, as in the general case one would have to compute the inverse function by means of zero-finding routines, which are typically slow and do not ensure convergence, especially when working in the complex plane. Consider also that the DL anomalous dimension from Eq. (4.1) is itself obtained by means of zero-finding routines, applied to the DL kernel which is given in terms of χ s , giving rise to nested zero-finding which clearly cannot guarantee best performance. Therefore, using the analytic expression Eq. (B.19) for χ s allows to have a single layer of zero-finding routines, improving significantly the numerical stability and the speed of the code. 20 The expansion of χ s in power of α s is given by Up to this order, these are identical to the dual of the exact anomalous dimension, as an obvious consequence to the fact that the approximate anomalous dimension is constructed to preserve LL and NLL accuracy.

B.4 Inverse Mellin transforms
Here we compute the various ingredients for the inverse Mellin transforms of Eqs. Because the HELL-x code, where these expressions are implemented, has to be fast, as it is meant to be used in PDF fits, the appearance of polylogarithms is not ideal. Therefore, we can consider a small-x approximation of these expressions. After all, the complicated structure of the O(α 3 s ) contribution in Eq. (4.31c) comes from the complicated all-order structure of the resummed result, but what really matters in the resummation is the prediction of small-x contributions, while uncontrolled terms which vanish as x → 0 are irrelevant. Hence, we approximate the expressions above as In these equations we have also provided the Mellin transform of the approximate expressions, which is needed for the analytic computation of the momentum conservation, Eq. (4.42). We verified that using these approximate expressions leads only to tiny deviations with respect to the exact expressions, and all in a region of x which is not under control of small-x resummation (specifically x > 10 −2 ). On the other hand, the speed-up is significant, fully justifying their use.