Towards parton distribution functions with small-$x$ resummation: HELL 2.0
Marco Bonvini Simone Marzani Claudio Muselli

=1 
             

 
                     

 
 










































































    ()
        [
    ]
    
    a]Marco Bonvini,[a]Dipartimento di Fisica, Sapienza Università di Roma and INFN, Sezione di Roma 1,
     Piazzale Aldo Moro 5, 00185 Roma, Italy
    b]Simone Marzani[b]Dipartimento di Fisica, Università di Genova and INFN, Sezione di Genova,
     Via Dodecaneso 33, 16146, Italy
    c]and Claudio Muselli[c]Tif Lab, Dipartimento di Fisica, Università di Milano and INFN, Sezione di Milano,
     Via Celoria 16, 20133 Milano, Italy
    https://arxiv.org/abs/1708.075101708.07510
    marco.bonvini@roma1.infn.itsimone.marzani@ge.infn.itclaudio.muselli@mi.infn.it
    
    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 , version ,
    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.
    
    
     Towards  parton distribution functions with small-x resummation: HELL 2.0
        [
        August 24, 2017
    ==========================================================================
    
    
    
    § 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^3LO <cit.>.
    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 <cit.>
    or various differential distributions in top pair production <cit.>
    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 <cit.>,
    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^3LL for a few observables <cit.>. 
    
    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 <cit.> (see also Ref. <cit.>) and performed dedicated studies that included threshold resummation in both coefficient functions and PDFs  in the context of the production of heavy supersymmetric particles <cit.>.
    
    It has to be noted that the inclusion of threshold resummation in PDF fits did not pose particular challenges because in the widely used  scheme the splitting functions that govern DGLAP evolution are not enhanced at large x <cit.>, 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.
    
    High-energy resummation of PDF evolution is based on the BFKL equation <cit.>.
    However, the correct inclusion of LL and NLL corrections to DGLAP splitting functions is far from trivial. 
    
    This problem received great attention in 1990s and early 2000s by various groups, see Refs. <cit.>, Refs. <cit.>, and Refs. <cit.> (for recent work in the context of effective theories, see <cit.>).
    
    High-energy resummation of partonic cross sections is based on the so-called -factorization theorem <cit.>, which has been used to compute high-energy cross sections for various processes:
    deep-inelastic scattering (DIS) <cit.>,
    heavy quark production <cit.>,
    direct photon production <cit.>,
    Drell-Yan <cit.>,
    and Higgs production <cit.>.
    The formalism has been subsequently extended to rapidity <cit.> and transverse momentum distributions <cit.>.
    
    Despite the wealth of calculations, it has proven very difficult to perform resummed phenomenology. Recently, some of us overcame these difficulties and developed a framework and a public code named  (High-Energy Large Logarithms) <cit.>, which is based on the formalism developed by Altarelli, Ball and
    Forte (ABF) <cit.>, but does contain significant improvements. 
    
    
    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 <cit.> 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 <cit.>,
    which provides coefficient functions with power-behaving 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  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 <cit.> that we inherited from the original ABF work <cit.>, 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  expansion.
    
    
    All these improvements are included in a new release of , version , which is publicly available at the webpage
    https://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^3LO. 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. <cit.>.
    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. <cit.>.
    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 <cit.>.
    
    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. <cit.>.
    
    
    
    This paper is organized as follows. In Sect. <ref> 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. <ref> we present our new realization of the resummation of running coupling contributions,
    and later in Sect. <ref> 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. <ref>, before concluding in Sect. <ref>. In Appendix <ref> we provide analytical results for the off-shell DIS partonic cross section with mass dependence,
    while in Appendix <ref> we collect some technical details on the actual implementation of the resummation in .
    
    
    
    
    
    § 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,
    
    F_a(x,Q^2)= x ∑_i∫_x^1  dz/z  C_a,ix/z, m_c^2/Q^2, m_b^2/Q^2,m_t^2/Q^2, f_i (z,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. <ref>.
    
    
    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=^2=Q^2, so =(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 ∼^2, 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 <cit.>;
    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. (<ref>) we consider Mellin moments of the structure functions
    
    F_a(N,Q^2) ≡∫_0^1 dx  x^N-1 F_a(x,Q^2)= ∑_i C_a,iN, m_c^2/Q^2, m_b^2/Q^2,m_t^2/Q^2, f_i(N,Q^2),
    
    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:
    
        C_a,iN, m_c^2/Q^2, m_b^2/Q^2,m_t^2/Q^2,   = ∫_0^1 d z   z^N  C_a,iz, m_c^2/Q^2, m_b^2/Q^2,m_t^2/Q^2,, 
        
        f_iN, Q^2   = ∫_0^1 d z   z^N  f_iz, Q^2.
    
    The last equation implies that the DGLAP anomalous dimensions are defined as
    γ_ij(N,)=  ∫_0^1 d z   z^N  P_ijz ,,
    
    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 ^n P^(n-1)_ij∼^n 1/zln^n-11/z. In Mellin space these logarithms are mapped into poles in N=0, which results in the following LL behaviour for the anomalous dimension: ^n γ^(n-1)_ij∼/N^n. In practice, not every entry of the anomalous dimension matrix is LL at small-x.
    
    On the other hand, the behaviour of the DIS coefficient functions in Mellin space is ^n C^(n)_a,i∼/N^n-1, 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 (^0) 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 <cit.> 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  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
    
    F_a= F_a^ S+ F_a^ NS.
    
    The singlet structure functions contain both gluon and quark (singlet) contributions
    
    F_a^ S = C_a,g  f_g + C^S_a,q  f_S ,
    
    where 
    
    f_S = ∑_i=1^n_f
        f_q_i+f_q̅_i
    ,
    n_f denoting the number of active quark flavours.
    The non-singlet structure function instead reads, see e.g. <cit.>,
    
    F_a^ NS = C^ NS_a,q ∑_i=1^n_fe_i^2f_q_i +f_q̅_i-1/n_ff_S,
    
    where e_i is the electric charge of the quark q_i, i.e. its coupling to the photon.[Here we assume that the DIS interaction is mediated by a photon, thus ignoring a possible
    contribution from the Z boson, or the charged-current case in which a W is exchanged.
    This choice makes the discussion here and in the following somewhat simpler, and allows to focus on the details of the resummation.
    Most of what follows does not depend on this assumption, as small-x resummation affects just the singlet:
    in particular, the small-x logarithmic content of the singlet and pure-singlet coefficients, Eq. (<ref>), is identical (for a discussion about small-x enhancements in the non-singlet sector see, for instance, Ref. <cit.>).
    
    The generalization to a generic vector boson is rather straightforward, and will be presented in Sect. <ref>.
    ]
    Furthermore, we can collect the terms proportional to the singlet PDF, obtaining
    
        F_a    = C_a,g  f_g + C^ PS_a,q  f_S
        + C^ NS_a,q ∑_i=1^n_fe_i^2f_q_i+f_q̅_i,
    
    where we have defined the so-called pure-singlet coefficient function,
    
    C^ PS_a,q = C_a,q^S - ⟨ e^2 ⟩ C^ NS_a,q,
    
    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- ∑_k=1^n ^k   C^res,(k),
    
    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 <cit.>.
    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. <cit.>.
    
    
    
    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 <cit.>.
    
    The PDFs in the two schemes are related by matching conditions
    
    f_i^ = ∑_j=g,q_1,q̅_1,…, q_n_f,q̅_n_f K_ij^  f_j^,     i=g,q_1,q̅_1,…, q_n_f,q̅_n_f, q_n_f+1,q̅_n_f+1,
    
    where the sum runs over active flavours in the n_f scheme, and K_ij^ are matching functions. 
    
    Note that we only consider here factorization schemes in which the matching coefficients depend on the heavy quark mass only through logarithms of Q^2/^2.
    In particular, we will consider only -like schemes, where all the PDFs f_i^ 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. (<ref>),
    have a purely perturbative origin, and depend effectively on the heavy quark mass.
    
    The structure functions Eq. (<ref>) can be written in either scheme
    
        F_a
             = F_a^≡ C^_a,g f^_g + C^ PS_a,q f^_S + C^ NS_a,q∑_i=1^n_fe_i^2f^_q_i+f^_q̅_i
           = F_a^≡ C^_a,g f^_g + C^ PS_a,q f^_S + C^ NS_a,q∑_i=1^n_f+1e_i^2f^_q_i+f^_q̅_i.
    
    To all orders in , 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. (<ref>):
    
        C^_a,g   = C^_a,g K_gg^ + C^ PS_a,q (n_fK_qg^+K_hg^) + C^ NS_a,q∑_i=1^n_fe_i^2 K_qg^ + e_n_f+1^2 K_hg^,
    
    and similarly for the quark coefficient functions.
    In the equation above, we have defined
    
    K_hg^≡ K_q_n_f+1g^ + K_q̅_n_f+1g^
    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. (<ref>) 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. <ref>.
    For a detailed and general discussion, not limited to small x, see e.g. <cit.>.
    
    
    
    
      §.§.§ 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. <cit.>, where explicit resummed results were presented in the DIS factorization scheme. 
    
    Here, we extend the results to -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. <ref>.
    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)
    Δ C_a,i^ = ∑_k=1^n_f e_k^2 _a,i(m_q_k) + ∑_k=n_f+1^6 e_k^2  _a,i(m_q_k),      i=g,q,
    
    where _a,i(m_q_k) is the resummed contribution in the case of a light active flavour being struck by the photon,
    and _a,i(m_q_k) is the resummed contribution in the case of a heavy flavour being struck by the photon.[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 _a,i and _a,i
    would also depend on the mass of the outgoing quark.]
    Recall that being active does not necessarily imply being massless and indeed both massless and massive flavours contribute to _a,i.
    Crucially though, in this contribution collinear logarithms are factorized and resummed and, consequently, the zero mass limit of _a,i is finite. 
    On the other hand, _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
    ∑_k=1^n_f e_k^2 _a,i(0) = ⟨ e^2 ⟩ n_f _a,i(0) ,      i=g,q.
    
    For this reason, a factor n_f is usually included in the definition of the massless singlet coefficient function,
    see e.g. Refs. <cit.>.
    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 <cit.> or FONLL <cit.>
    (which are formally identical <cit.>)
    in which the mass dependence is retained in the coefficient functions.
    
    
    
    We now perform a logarithmic counting on Eq. (<ref>).
    None of the matching functions are LL, with the exception of the LO diagonal components, which are all equal to 1 at (^0), K_ii^(0)=1.
    Furthermore, all coefficient functions are NLL, except the non-singlet LO coefficient of F_2,
    which is C^ NS,(0)_2,q=(1) and thus LL.
    The leading non-trivial logarithmic contributions in the coefficient functions are then NLL,
    and their resummed contributions, Eq. (<ref>), are related in the two schemes by
    
        Δ C^_a,g   = Δ C^_a,g + e_n_f+1^2 C^ NS,(0)_a,qΔ K_hg(), 
        Δ C^_a,q   = Δ C^_a,q + e_n_f+1^2 C^ NS,(0)_a,qΔ K_hq(),
    
    where Δ K_ij are the NLL resummed contributions to K_ij^.
    In the results above, we have neglected all contributions which are products of two NLL functions.
    We have also dropped K_qg^ (and K_qq^), 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 .
    Note that, for simplicity, we are not indicating in Δ K_ij the label ^, but we emphasize its (logarithmic) dependence on the heavy quark mass.
    
    Note that C^ NS,(0)_a,q in Eq. (<ref>) 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. (<ref>), is undetermined: there are two more unknowns than equations <cit.>.
    This choice is the one leading to S-ACOT/FONLL <cit.>
    and TR <cit.> schemes, where we have
    
        C^ NS,(0)_2,q(0)   =1, 
        
        C^ NS,(0)_L,q(0)   =0.
    
    In alternative approaches, like ACOT <cit.> or, equivalently, a new incarnation of FONLL <cit.>, 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:
    
        C^ NS,(0)_2,q()
           =√(1+4^2/Q^2)2/1+√(1+4^2/Q^2)^N+1, 
        
        C^ NS,(0)_L,q()
           =4^2/Q^2/√(1+4^2/Q^2)2/1+√(1+4^2/Q^2)^N+1.
    
    Note that the massless limit of these expressions reduces to the massless result, Eq. (<ref>).
    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. <ref>.
    
    We now consider again the decomposition Eq. (<ref>).
    In the n_f+1 scheme it simply becomes
    Δ C^_a,i = ∑_k=1^n_f+1 e_k^2 _a,i(m_q_k) + ∑_k=n_f+2^6 e_k^2  _a,i(m_q_k),
         i=g,q,
    
    where _a,i() 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. (<ref>) and (<ref>) into Eq. (<ref>).
    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:
    
    e_n_f+1^2  _a,i() = e_n_f+1^2  _a,i() + e_n_f+1^2 C^ NS,(0)_a,qΔ K_hi(),
         i=g,q.
    
    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. <ref>.
    We then find the (collinearly factorized) massive coefficients in the n_f+1 scheme for F_2,
    assuming Eq. (<ref>),
    
        _2,g()    = _2,g() - Δ K_hg(), 
        _2,q()    = _2,q() - Δ K_hq(),
    
    and for F_L
        _L,g()    = _L,g(), 
        _L,q()    = _L,q().
    
    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. (<ref>) 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^_ij (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. <cit.>), 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 _a,i(),
    such that the massless limit →0 of the coefficient functions in the new scheme _a,i() is finite.
    If we further require that the n_f and n_f+1 scheme are both of the same type, e.g. -like,
    we also need to impose that the massless limit of the massive coefficient _a,i()
    is just what we would have computed if the (n_f+1)-th flavour were massless:
    lim_Q≫_a,i() = _a,i(0),      i=g,q.
    
    (Note that, at the logarithmic accuracy we are interested in,
    _a,i(0) are the same in the n_f and in the n_f+1 schemes.[In the actual resummed expressions we compute, there is a non-zero n_f dependence in the coefficient functions _a,i(0),
    which is however subleading and can therefore be ignored.])
    The massless limit Eq. (<ref>) 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 _a,i(m_q_n_f+1).
    It does not tell us anything about the power corrections in /Q in the n_f+1 scheme, which have to be determined
    by the matching procedure.
    
    The results Eq. (<ref>) 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. (<ref>), 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≫.
    
    We can exploit the last consideration to derive the desired resummed expressions for the matching functions Δ K_hi, i=g,q.
    Indeed, Eq. (<ref>) can be inverted to give Δ K_hi in terms of the massive coefficient functions _2,i(m),
    which are known, and of the coefficient functions _2,i(m) in the n_f+1 scheme, the massless limits of which, _2,i(0), are known.
    We can therefore consider the massless limit →0 of this expression, keeping in mind that Δ K_hi and _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  or mass independent terms, we then find
    Δ K_hi() = lim_Q≫_2,i() - _2,i(0) ,      i=g,q,
    
    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 _2,i(m) through Eq. (<ref>),
    as well as the PDFs in the n_f+1 scheme according to Eq. (<ref>).
    However, this “factorization” of power behaved contributions is beyond the control of the collinear factorization framework,
    so the simplest choice Eq. (<ref>) is perfectly acceptable and therefore universally adopted.
    
    The definition Eq. (<ref>) 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[These colour-charge relations are strictly speaking only valid when the resummed contributions refer to the pure-singlet.
    If they refer to the singlet, in the massless case there is a fixed-order contribution which needs to be subtracted <cit.>.
    Alternatively, we can say that these relations are valid for k_a,i with k≥1.]_a,q(0) = C_F/C_A_a,g(0),     _a,q(m) = C_F/C_A_a,g(m).
    
    From Eq. (<ref>) it then immediately follows
    Δ K_hq(m) = C_F/C_AΔ K_hg(m).
    
    Together, Eqs. (<ref>) and (<ref>) imply through Eqs. (<ref>) and (<ref>)
    a colour-charge relation
    _a,q(m) = C_F/C_A_a,g(m).
    
    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_V^2+g_A^2, 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_V^2+g_A^2, but a contribution
    proportional to g_A^2/(g_V^2+g_A^2) 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. <ref>, remains unchanged, provided the photon couplings
    e_k^2 are replaced with the more general coupling (g_V^2+g_A^2)_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-quark q̅',
    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. (<ref>) is different for q and q̅, depending on the process.
    In particular, this implies that when using Eq. (<ref>) to derive Eq. (<ref>),
    only either the heavy quark q_n_f+1 or the heavy anti-quark q̅_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^ CC, NS,(0)_3,q=1 and C^ CC, NS,(0)_3,q̅=-1.[Note that the opposite sign in these two contribution makes the collinear subtractions
    in the fully massless case ineffective. Indeed, the collinear singularities in this case cancel automatically.
    In fact, after cancellation, the non-singular part vanishes, thus making the massless singlet contribution to F_3 zero
    (which remains true at higher orders, since the underlying reason is the antisymmetry of the F_3 contribution for the exchange of the two quark masses).
    The same mechanism holds in the neutral-current massive case.]
    Putting everything together, we have that the analogous of Eqs. (<ref>) and (<ref>) are
    
        _2,i^ CC()    = _2,i^ CC() - Δ K_hi()/2,          i=g,q, 
        _L,i^ CC()    = _L,i^ CC(), 
        _3,i^ CC()    = 
                                        _3,i^ CC() - Δ K_hi()/2    if final state is q Q̅,
        _3,i^ CC() + Δ K_hi()/2    if final state is Q q̅,
    
    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, _3,i^ CC(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. <ref> 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. (<ref>), we also determine the all-order behaviour of the matching functions.
    
    The all-order behaviour of partonic coefficient functions is obtained using the -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. <cit.>. 
    
    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 <ref>, where we collect all the relevant results for the off-shell DIS coefficient functions.
    These results have been implemented in the public code , version ,
    and are also accessible through the public code  <cit.>, to which  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. <cit.>.
    Following Ref. <cit.>, in our recent work <cit.> 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. <ref> for details about our notation),
    is convoluted with an evolution operator U(N,ξ),
    
        0_2,g(0)    = -∫ dξ d/dξ C_2(0,ξ,0,0)  U(N,ξ) - U_qg/n_f, 
        0_L,g(0)    = -∫ dξ d/dξ C_L(0,ξ,0,0)  U(N,ξ),
    
    where
    
        U(N,ξ)= exp
            ∫_1^ξd ζ/ζγ_+ N, (ζ Q^2)
        ,
    
    and γ_+ is the small-x resummed anomalous dimension.
    
    Note that the above expressions hold in a factorization scheme denoted Q_0 <cit.>, which differs from  at relative (^3). In the context of small-x resummation this scheme is preferred because it leads to more stable results <cit.>.
    
    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 <cit.>
        U_qg = ∫dξ/ξγ_qg(N,(ξ Q^2)) θ(1-ξ)   U(N,ξ)
    
    where γ_qg is the resummed qg anomalous dimension.
    All ξ integrals extend to ∞ and start from the position of the Landau pole, ξ_0=exp-1/β_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. <cit.>).
    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 <ref> 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,
    0_a,g(m) = -∫ dξ d/dξ C_a0,ξ,ξ_m,ξ_m  U(N,ξ),    a=2, L,
    
    where we have defined ξ_m=m^2/Q^2 (see again App. <ref> 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 <cit.>
    (see also Ref. <cit.>).
    We report them in Appendix <ref>, 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. <ref> 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
    
        0_2,g^ CC(m)    = -∫ dξ d/dξ C_20,ξ,ξ_m,0  U(N,ξ) - U_qg/2n_f, 
        0_L,g^ CC(m)    = -∫ dξ d/dξ C_L0,ξ,ξ_m,0  U(N,ξ) - ξ_m/1+ξ_m U_qg/2n_f, 
        0_3,g^ CC(m)    = 
           -∫ dξ d/dξ C_30,ξ,ξ_m,0  U(N,ξ) + 1/1+ξ_m U_qg/2n_f    (qQ̅), 
         -∫ dξ d/dξ C_30,ξ,0,ξ_m  U(N,ξ) - 1/1+ξ_m U_qg/2n_f    (Qq̅).
    
    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. (<ref>).
    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. <ref>.
    Note that, according to the notation defined in App. <ref> 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 off-shell 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 <ref>.
    
    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. (<ref>) (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. (<ref>), (<ref>),
    for F_2 and F_3 one further needs to subtract the matching condition, Eqs. (<ref>), (<ref>), (<ref>).
    Accordingly, the resummed massive coefficient functions for the massive active flavour are given in neutral current by
    
        0_2,g(m)    = -∫ dξ d/dξ C_20,ξ,ξ_m,ξ_m  U(N,ξ) - Δ_0 K_hg(m),
        0_L,g(m)    = -∫ dξ d/dξ C_L0,ξ,ξ_m,ξ_m  U(N,ξ),
    
    and in charged current by
    
        0_2,g^ CC(m)    = -∫ dξ d/dξ C_20,ξ,ξ_m,0  U(N,ξ) - U_qg/2n_f - Δ_0 K_hg(m)/2,
        0_L,g^ CC(m)    = -∫ dξ d/dξ C_L0,ξ,ξ_m,0  U(N,ξ) - ξ_m/1+ξ_mU_qg/2n_f, 
        0_3,g^ CC(m)    = 
           -∫ dξ d/dξ C_30,ξ,ξ_m,0  U(N,ξ) + 1/1+ξ_mU_qg/2n_f - Δ_0 K_hg(m)/2    (qQ̅), 
         -∫ dξ d/dξ C_30,ξ,0,ξ_m  U(N,ξ) - 1/1+ξ_mU_qg/2n_f + Δ_0 K_hg(m)/2    (Qq̅)
    
    (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. (<ref>) (and equivalently Eq. (<ref>)) with Eq. (<ref>) 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_a,g^ (CC)(m), a=2,3 does not commute with the ξ integration.
    
    
    
    
      §.§.§ Resummed matching functions
    
    
    We can now use Eq. (<ref>) to determine the resummed matching function Δ_0 K_hg(m).
    Using Eqs. (<ref>) and (<ref>) we can write[Note that this equation can be written in two alternative forms by comparing Eq. (<ref>) to Eq. (<ref>)
    or by noting that the massless limit of Eq. (<ref>) vanishes.
    Using the results of App. <ref> it is easy to verify that all these forms are equivalent and lead to the same result.]
        Δ_0 K_hg(m)    = U_qg/n_f
        +∫ dξ d/dξ C_2(0,ξ,0,0)  U(N,ξ)
          -lim_ξ_m→0∫ dξ d/dξ C_20,ξ,ξ_m,ξ_m  U(N,ξ),
    
    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/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 <ref>)
    Δ_0 K_hg^f.c.(m) = U_qg/n_f-
    /πm^2/Q^2^γ_s1-γ_s/γ_sΓ^3(1-γ_s)Γ(1+γ_s)/(3-2γ_s)Γ(2-2γ_s).
    
    A non-trivial check of the above expression can be performed by considering its perturbative expansion,
    
        Δ_0 K_hg^f.c. (m)
             = [h_0/n_f/γ_s+h_1/n_f+(h_2/n_f)γ_s+…
                  -1/3πγ_s - 5+3log(m^2/Q^2)/9π - 56+30log(m^2/Q^2)+9log^2(m^2/Q^2)/54πγ_s +…] 
           = - log(m^2/Q^2)/3/π- 28+30log(m^2/Q^2)+9log^2(m^2/Q^2)/18N/π^2 +(^3),
    
    where h_i are the (known) coefficients of the perturbative expansion of γ_s U_qgf.c.=γ_qg
    in powers of γ_s <cit.>.
    The first coefficients read h_0=n_f/3π, h_1=n_f/3π5/3, h_2=n_f/3π14/9.
    Note that the collinear pole 1/γ_s cancels out in the sum.
    The second equality is then obtained replacing γ_s= C_A/π N+(^2).
    
    The above expression can be compared with the fixed-order results presented in Ref. <cit.>.
    The  term corresponds to the Mellin transform of the NLO result computed in N=0,
    and the ^2 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 to , which affects the second term in Eq. (<ref>), we find full agreement with the high-energy limit of the ^3 result <cit.>.
    
    We can now restore the running-coupling effects in the resummation from the fixed-coupling result by computing
    Δ_0 K_hg(m) = U_qg/n_f - ∫ dξ d/dξ K_hg(ξ,ξ_m)  U(N,ξ),
    
    where 𝒦_hg is obtained as the inverse Mellin transform of its fixed-coupling counterpart, second term in Eq. (<ref>).
    The computation of this inverse Mellin transform is done in App. <ref>, and its derivative is given by
    
        d/dξ𝒦_hg(ξ,ξ_m)
           = /3π 6ξ_m/ξ^2
            1 - 4ξ_m/ξ√(ξ/ξ+4ξ_m)log√(ξ/4ξ_m)+√(1+ξ/4ξ_m)
        .
    
    Clearly, by construction, Eq. (<ref>) with Eq. (<ref>) reproduces the correct result
    in the fixed-coupling limit, Eq. (<ref>).
    Also, it clearly includes the correct resummation of the subleading running coupling effects, as the form Eq. (<ref>)
    is the standard expression for such resummation <cit.>.
    Eq. (<ref>) is a new result. It allows to resum the matching conditions in -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 () and (^2) contributions to the resummed matching function are explicitly given in Eq. (<ref>)
    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 U_qg/n_f, which is in turn needed for the massless collinear subtractions
    of DIS coefficient functions in Eqs. (<ref>), (<ref>) and (<ref>).
    To complete the list, one needs to expand the ξ-integrals of the off-shell coefficient functions with the evolution factor U(N,ξ).
    As for the matching function, up to (^2) this expansion can be simply obtained by working in the fixed-coupling limit.
    For this, we need the expansions in M of the Mellin transforms of the off-shell coefficients, given in App. <ref>,
    where M should be replaced by γ_s= C_A/π N+(^2) and expanded out.
    By doing so, all the 1's and 2's can be constructed, making   the matching of each ingredient with the corresponding fixed order up to NNLO straightforward.
    
    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. <ref>.
    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 <cit.>)
    S-ACOT construction, with a canonical perturbative counting based on explicit powers of 
    (at fixed /N for resummed contributions).
    This is equivalent to a plain (i.e., without damping <cit.>) 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  package,
    with which our results can be directly used for resummed DIS phenomenology.
    
    The FONLL approach <cit.> 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 <cit.> of structure functions, assuming a single heavy quark q_n_f+1 with mass m, is performed as
    
    F_a^ FONLL(m) = F_a^(m) + F_a^(0) - F_a^ d.c.(m),      a=2,L,3,
    
    where F_a^(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_a^(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_a^ d.c.(m) is the double-counting term (called massive-zero),
    which can be either seen as the fixed-order expansion of F_a^(0) or the “massless limit”
    (in which divergent terms are kept finite) of F_a^(m).
    The combination F_a^(0) - F_a^ d.c.(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) - F_a^ d.c.(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.
    F_a^ FONLL(m) = F_a^(m).
    
    Thus, the double counting term, which is the only non-trivial ingredient in Eq. (<ref>),
    is given by
    
    F_a^ d.c.(m) = F_a^(m) + F_a^(0) - F_a^(m)
    
    The corresponding small-x resummed coefficient functions, analogous to _a,i and _a,i,
    can be obtained from Eq. (<ref>) using Eq. (<ref>)
    together with Eqs. (<ref>) and (<ref>),
    and are given by
    ^ d.c._a,i(m) = _a,i(0) + C^ NS,(0)_a,qΔ K_hi(m),      i=g,q
    
    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. (<ref>),
    we can also consider the version of FONLL which includes a damping on the resummed contribution,
    
    F_a^ FONLL+damp(m) = F_a^(m) + θ(1-ξ_m)1-ξ_m^2
        F_a^(0) - F_a^ d.c.(m)
    ,
    
    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 _a,i(m).
    
    A word of caution is needed when discussing the perturbative counting. The canonical counting would consist in including
    all contributions at () for NLO and all contributions at (^2) 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_a^(0) is retained at (), as well as the matching functions,
    but the massive contribution F_a^(m) is computed at one order higher, (^2) <cit.>.
    
    When this particular perturbative counting is adopted, the double-counting piece must be computed with care.
    In particular, only the definition of F_a^ d.c.(m) as the fixed-order expansion of F_a^(0)
    to (^2) gives the correct result.
    
    As far as small-x resummation is concerned, one has to use 2_a,i(m) for the massive part
    and 1_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_a,i(0) and Δ_1 K_hi(m) in Eq. (<ref>) is needed.
    It can be explicitly verified that in this way the “resummed contribution”F_a^(0)-F_a^ d.c.(m)
    is indeed of (^3) 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 <cit.> construction, and it
    also corresponds to the variant FONLL_IC proposed in Refs. <cit.>
    to account for a possible intrinsic component of the charm PDF.
    Following the latter references, we define
    
        δ F_a
           = F_a^ FONLL_ IC - F_a^ FONLL
           = F_a^ ACOT - F_a^S-ACOT,
    
    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. (<ref>),
    and the one obtained with massless NS coefficients, Eq. (<ref>).
    We thus have (we apologize with the Reader for the awkward notation)
    
        Δδ c_a,i(m)    = 
            C^ NS,(0)_a,q(0) - C^ NS,(0)_a,q(m)
        Δ K_hi(m),       a=2,L,       i=g,q, 
        Δδ c_a,i^ CC(m)    = 
            C^ CC,NS,(0)_a,q(0) - C^ CC,NS,(0)_a,q(m)
        Δ K_hi(m)/2,       a=2,L,3,       i=g,q,
    
    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. (<ref>), 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 <cit.>,
    which we followed with minor modifications in our previous work <cit.>,
    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 <cit.> 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. <cit.>),
    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 -dependence and the non-trivial M-dependence of the DL kernel.
    In some approaches, e.g. <cit.>,
    the equation is thus solved numerically, and the resummed anomalous dimension derived in a numerical way.
    Instead, in Refs. <cit.> an approximate analytic solution,
    in which both - 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. <cit.>
    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  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
    χ_ DLγ_ DL(N,), = N,
    
    where χ_ DLM, is the BFKL kernel and γ_ DL(N,) the DGLAP anomalous dimension.
    The knowledge of the BFKL kernel at N^kLO provides by duality all the N^kLL 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^kLO,
    thus implying that they also contain (by duality with each other) all N^kLL contributions.
    Therefore, dual DL kernels obtained with fixed N^kLO in both, usually denoted DL-N^kLO, are matched results
    of the form N^kLO+N^kLL, and so they are both double (next-to-^k) 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. (<ref>) assumes that the strong coupling  does not run, namely it is Q-independent.
    When the running of  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 /N
    times an overall factor  <cit.>γ_ss^ rc(N,) = -β_0.χ_0”(M)χ_0(M)/2χ_0'^2(M)|_M=γ_s(/N),
    
    where γ_s(/N) is the dual of the LO BFKL kernel χ_0(M).
    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(/N)=1/2, i.e. N=χ_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 /N times a function of β_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  is not fixed but is running.
    When the coupling runs, (Q^2) becomes a differential operator =/(1-β_0 d/dM) in Mellin space,
    with =(Q_0^2), Q_0 being a fixed scale, and where we have assumed 1-loop running.
    In principle, there can be  computed at different scales in the kernel, but one can always rewrite it as
    (Q^2) by evolving to that scale and expanding the relevant evolution factors.
    In this way, the  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  can be set equal to (Q^2) in all terms without
      modifying explicitly the kernel, at DL-NLO changing the argument of each  to (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. (<ref>) later)
      does not correspond, by construction, to a kernel in which all 's are (Q^2).
      Therefore, at NLO, the form of the kernel in which all 's are (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. <cit.>.
    
      *  Once all powers of  are computed at Q^2,
      which is equivalent to say that all powers of  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  is retained
      and all others are frozen to (Q_0^2), is used.
    
    The fact that the complicated all-order  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 β_0 times a LL function of /N,
    which originates from 1-loop running at lowest order. The NLO contribution to the BFKL kernel would produce corrections
    which are of order (β_0)^n (/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. <cit.> and used later in Ref. <cit.>
    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 <cit.>,
    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. <ref>.
    
    Having understood that the linear  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. <cit.>, 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 -expansion of the resulting anomalous dimension contains
    half-integer powers of .
    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 . 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  approximation and collinear M approximation.
    The starting point is the on-shell BFKL kernel in symmetric variables <cit.> that we denote simply χ(M,),
    whose  dependence is approximated as
    χ(M,)
    = χ(M,) + (-) χ'(M,)
    = χ̅(M,) + χ'(M,)
    
    where prime denotes derivative with respect to .
    It is important to observe that this approximation includes an (^0) term,
    χ̅, which is not physical and not present in the original kernel which is of lowest order ().
    For this reason the kernel Eq. (<ref>) does not go to zero as →0.
    One could therefore consider another linear approximation,
    χ(M,)
    = χ(M,)/,
    
    which does go to zero as →0, but does not reproduce the exact derivative in =.
    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,)=χ_0(M).
    In the following, we will use the first approximation, Eq. (<ref>), to derive our results,
    which can then be easily translated into the second approximation, Eq. (<ref>), by simply letting χ̅=0, χ' = χ/.
    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. (<ref>),
    from which the resummed anomalous dimension can be derived,
    is given by <cit.>
        N - χ̅(M,)
     f(N,M) = χ'(M,) f(N,M),
    
    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
    γ_ rc(N,) = 
        d/dtlog∫_M_1-i∞^M_1+i∞dM/2π i e^Mtexp∫_M_0^M dM' 1/β_01-χ'(M',)/N - χ̅(M',)
    _t=0,
    
    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,
    χ(M,) = c() + κ()/2M-M_ min()^2  + (M-M_ min)^3,
    
    where M_ min differs in general from 1/2 by terms of (). 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  in momentum space.
    A better approximation, which we propose here, is inspired by a collinear plus anti-collinear approximation of the kernel <cit.>
    generalized to account for a minimum which is not in 1/2:
    χ_ coll(M,) = A()
        1/M+1/2M_ min()-M
     + B().
    
    Expanding around its minimum M=M_ min we find
    χ_ coll(M,) = B+2A/M_ min + 2A/M_ min^3M-M_ min^2 + (M-M_ min)^3,
    
    which leads to the identifications
    
    A=M_ min^3κ/4,      B=c - M_ min^2κ/2,
    
    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  dependence of the collinear kernel Eq. (<ref>).
    For A() and B(), we use the very same linear decomposition Eq. (<ref>),
    while the position of the minimum M_ min() will not be considered as an operator.[Note that this assumption is less crude than the approach of previous works, where M_ min was simply approximated to be 1/2 in the RC equation.]
    The integrals in Eq. (<ref>) 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. (<ref>) is given by
    χ'(M,)/N - χ̅(M,)
    = c'+κ'/2(M-M_ min)^2/N-c̅-κ̅/2(M-M_ min)^2,
    
    while for collinear kernel we find
    
        χ'(M,)/N - χ̅(M,)   = 2M_ minA'+B'M(2M_ min-M)/(N-B̅)M(2M_ min-M)-2M_ minA̅
           = c'-B'/M_ min^2(M-M_ min)^2/N -c̅ - (N-B̅)/M_ min^2(M-M_ min)^2,
    
    having used in the last step Eq. (<ref>). By direct comparison, we find the translation rules
    κ'→ -2B'/M_ min^2 = κ'-2c'/M_ min^2,    κ̅→ 2N-B̅/M_ min^2 = κ̅+2N-c̅/M_ min^2.
    
    Hence, the final solution is given by <cit.>[In order to account for a generalized position of the minimum, we have recomputed the solution analytically,
    thus providing a useful cross-check of the result presented in the literature.]γ_ rc(N,) = M_ min + β_0
        z k_ν'(z)/k_ν (z) -1
    ,
    
    where k_ν(z) is a Bateman function, with
    
        1/   = 1/ + κ'-2c'/M_ min^2/κ̅+2(N-c̅)/M_ min^2
        
          z    = 1/β_0√(N-c̅/κ̅/2 +(N-c̅)/M_ min^2)
        ν   = c'/N-c̅ + κ'-2c'/M_ min^2/κ̅+2(N-c̅)/M_ min^2 z.
    
    We immediately observe that the limit c̅,κ̅→0 of these expressions is finite and trivial,
    so the solution in the approximation Eq. (<ref>) 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. (<ref>) with Eq. (<ref>) represents a new result with respect to the
    old “Bateman” RC solution of Refs. <cit.> 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. (<ref>),
    which is equivalent to an expansion in powers of β_0 of Eq. (<ref>).
    This expansion is also needed to identify the proper double counting with the DL part.
    We find
    
        γ_ rc(N,)    = M_ min - √(N-c/κ/2 +(N-c)/M_ min^2) - β_0
        +1/4β_0^2
            3κ'-2c'/M_ min^2/κ+2(N-c)/M_ min^2-c'/N-c
        
           +β_0^2^3 √(κ/2 +(N-c)/M_ min^2/N-c)c'κ-κ'(N-c)/32(N-c)^2[2(N-c)/M_ min^2+κ]^2
               ×[16c^2/M_ min^2 + 8N(κ+2N/M_ min^2) - c(16 c' /M_ min^2+8κ-3κ'+32N /M_ min^2) 
                    +(5c'κ+16c'N /M_ min^2-3κ'N)]+ (β_0^3).
    
    This result shows a number of interesting features, especially when compared with the analogous expansion in the case of the quadratic kernel <cit.>.
    First, we note that at large N all terms go to zero except for the running coupling correction -β_0, 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_ min^2 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 (^0), while it would be of ()
    if the denominator were just κ/2, as it happens with the quadratic kernel. Having a denominator of ()
    produces an  expansion of this result which contains half-integer powers of .
    
    Instead, the  expansion of Eq. (<ref>), and hence of Eq. (<ref>), is perfectly acceptable with only integer powers of .
    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. (<ref>),
    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  expansions of the RC contributions which will be needed in Sect. <ref> 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. (<ref>) have to be subtracted
    (the first two because they are LL, and the third because it is of (), and hence already included in the DL-LO).
    After subtraction we have the expansion
    
        Δ_DL-LOγ_ rc(N,)
           ≡γ_ rc(N,) - 
            M_ min - √(N-c/κ/2 +(N-c)/M_ min^2) - β_0
        
           = β_0^2 3κ_0/32-c_0/N + (^3),
    
    where κ_0 and c_0 are the derivatives of κ and c computed in =0, and are given by
    
    c_0 = C_A/π4log2,
        κ_0 = C_A/π28ζ_3.
    
    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,) <cit.>.
    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 (^2) of γ_ rc.
    However, we observe that at (β_0) the expansion of γ_ rc, Eq. (<ref>),
    contains terms which are formally NLL, and specifically given by β_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 <cit.>, Eq. (<ref>),
    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 <cit.>,
    
        Δγ_ss^ rc(N,)
           = -β_0
            χ_0”(M)χ_0(M)/2χ_0'^2(M) - 1
        _M=γ_s(/N)
           
        = (^4),
    
    where χ_0(M) is the LO BFKL kernel and γ_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 (^4). 
    Since the kernel used in the RC resummation is only approximate, the function γ_ rc does not correctly
    predict all the NLL contributions of Eq. (<ref>). Therefore, Eq. (<ref>) must be still added to the DL-NLO
    result, and the (β_0) part of γ_ rc
    has to be considered as a double counting term with respect to Δγ_ss^ rc, and hence subtracted.
    Thus, for RC resummation matched to DL-NLO, we further need to subtract the fourth term of Eq. (<ref>),
    
        Δ_DL-NLOγ_ rc(N,)
           ≡γ_ rc(N,) - [M_ min - √(N-c/κ/2 +(N-c)/M_ min^2) - β_0
                         +1/4β_0^23κ'-2c'/M_ min^2/κ+2(N-c)/M_ min^2-c'/N-c]
        
           = β_0^2^3 κ_0/16N + (^4),
    
    where we have also written its  expansion in the last line.
    We observe that, formally, only the (β_0) times a LL function is really doubly counted,
    so in principle one could expand the last line of Eq. (<ref>) 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. (<ref>) and (<ref>) at lowest order only involve lowest order derivatives of c and κ,
    Eq. (<ref>), which in turn are determined from just the LO BFKL kernel,
    and hence do not depend on the actual construction of the DL kernel.
    Also, they are the same irrespectively of the kind of approximate dependence on  is used,
    either Eq. (<ref>) or Eq. (<ref>).
    
    
    
    
    
    
     §.§ Singularity mismatch
    
    
    The functions Δ_DL-(N)LOγ_ rc(N,), Eqs. (<ref>) and (<ref>),
    contain the resummation of β_0 contributions which should be added to γ_DL-(N)LO, respectively.
    The square root term in the subtractions of Eqs. (<ref>) and (<ref>)
    contains the LL singularity which has to be removed from the DL result, and replaced with the pole singularity
    contained in the γ_ rc term.
    This cancellation is exact if the parameters of the minimum, c() and κ(), are computed
    from the same DL kernel used for the fixed-coupling duality which defines γ_DL-(N)LO.
    
    
    
    
    At LO+LL, the kernel is the one obtained putting on-shell Eq. (<ref>), so the singularity automatically cancels.[In fact, in Ref. <cit.> two different expressions of χ_s were used for computing the DL kernel
    and for the kernel used in RC resummation. Specifically, in the second case we used the dual of the exact LO anomalous dimension,
    which however could not be used for the DL kernel as the exact LO anomalous dimension γ_+ has a square-root branch-cut,
    due to the way the eigenvalue of the singlet anomalous dimension matrix is computed, which would produce a spurious
    oscillating behaviour. For this reason, for the DL we used an approximate LO anomalous dimension, thereby creating a mismatch
    in the singularities even at LO+LL.
    Here, thanks to the approximation discussed in App. <ref>, we use exactly the same kernel.]
    In Ref. <cit.> 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
    γ_ match^ LO+LL'(N,) = γ_ mN,c^ NLO(),κ^ NLO(),M_ min^ NLO()
    - γ_ mN,c^ LO(),κ^ LO(),1/2,
    
    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_ min^ (N)LO 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. (<ref>).
    
    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
    γ_ m(N,c,κ,M_ min) = M_ min - √(N-c/κ/2+(N-c)/M_ min^2) - M_ min^3κ/4N,
    
    has the desired properties. We note that, in contrast with the case of the quadratic kernel <cit.>, 
    here we do not need any further contribution, as this function already vanishes as 1/N at large N.
    Since the parameters in Eq. (<ref>) start differing at (^2), the function γ_ match^ LO+LL'
    is formally NLL. This is not a problem here, as the formal accuracy of the LO+LL^' result is just LL.
    The expansion of the matching function Eq. (<ref>) in powers of  is
    γ_ match^ LO+LL' (N,) = (^3).
    
    Hence, the choice of subleading terms in Eq. (<ref>) has the additional benefit that we do not need to keep
    the matching function into account when expanding the LO+LL^' result to (^2).
    
    
    At NLO+NLL, in Refs. <cit.> the kernel for RC resummation
    was constructed with a different  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. <ref>
    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 (^2) 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. <cit.>.
    We have verified that it is not possible to modify the function γ_ m, Eq. (<ref>),
    to make the function γ_ match NNLL without introducing new (uncanceled) singularities.
    To prove this statement, we consider a generalization of Eq. (<ref>)γ_ m(N,c,κ,M_ min) = M_ min - √(N-c/κ/2+(N-c)/M_ min^2) + η(N,c,κ,M_ min),
    
    where η(N,c,κ,M_ min) is a function to be determined,
    with the requirement that it must not introduce further leading singularities.
    Expanding Eq. (<ref>) to NLL, i.e. expanding in powers of  up to () at fixed /N, we find
    
        γ_ m(N,c,κ,M_ min)
             = 1/2 - √(N/-c_0/κ_0/2+4(N/-c_0)) + ηN/,c_0,κ_0,1/2
           +[m_1+c_1 κ_0 - c_0 κ_1 - 32 c_0^2 m_1 + κ_1 N/ + 64 c_0 m_1 N/ - 32 m_1 (N/)^2/√(2)√(N/-c_0)κ_0+8(N/-c_0)^3/2
               + c_1∂_2+ κ_1∂_3+ m_1∂_3ηN/,c_0,κ_0,1/2] + NNLL,
    
    where the index in the derivatives indicates with respect to which argument the derivative is computed,
    and m_1 is the () 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 (^2), 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= 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 Δγ_ss^ rc,
    Eq. (<ref>), and the Δ_DL-NLOγ_ rc, Eq. (<ref>).
    The latter exhibits explicitly a pole in N=c, which is different from an analogous singularity
    in the former,
    Δγ_ss^ rc(N,) ∼ -1/4β_0^2c_0/N- c_0,
    
    which is in N= 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. (<ref>) with Eq. (<ref>).
    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
    γ^ss_ match(N,) = 1/4β_0^2 
        c_0/N- c_0 - c'/N-c +c'-c_0/N
     = (^4),
    
    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 (^4) and therefore it does
    not contribute to the -expansion of the NLO+NLL result to (^3).
    Note that this singularity mismatch was present also in the original works using a quadratic kernel <cit.>;
    the problem there was solved by replacing by hand the singularity in Δ_DL-NLOγ_ rc
    with that of Δγ_ss^ rc, which effectively corresponds to using the same matching function Eq. (<ref>) but without the last term.
    
    
    
    
    
    
    
    
    
    § RESUMMED DGLAP EVOLUTION MATCHED TO NNLO
    
    
    As already discussed in Sect. <ref>, 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  expansion up to (^3), and replace it with the exact NNLO expression.
    While subtracting the NLO from the NLO+NLL is trivial, further subtracting the (^3) 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[Some developments with respect to our previous implementations have been performed, which make the code faster and more reliable.
      See App. <ref> for further details.]).
    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. <ref>.
    Indeed, for instance, this exercise allowed us to find a small mistake in the original ABF construction
    of the DL part <cit.>, which we also have
    inherited in our previous work <cit.>, and which we correct here.
    Then, in Sect. <ref>, 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 , version .
    
    
    
    
    
    
     §.§ Expansion of the Double Leading anomalous dimension
    
    
    In the ABF construction, the DL resummed anomalous dimension γ_ DL,
    Eq. (<ref>), is obtained from an implicit equation of the form
    χ_Σγ_ DL(N,),N, = N,
    
    where the function χ_ΣM,N, is a so-called off-shell BFKL kernel <cit.>.
    The DL anomalous dimension obtained through Eq. (<ref>) assumes fixed coupling ,
    and it thus receives a correction, Eq. (<ref>), due to running coupling effects.
    This correction starts at (^4) and it is therefore of no interest for the expansion of the DL result to (^3).
    In the following, we explain how to construct a perturbative expansion of γ_ DL(N,)
    defined by the implicit equation (<ref>), 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 (^2), which would be needed to match LL resummation to NLO.
    The off-shell kernel at LO, needed for LL resummation, is given by <cit.>χ_Σ^ LO(M,N,)
    = χ_s/M + χ_s/1-M+N + χ̃_0(M,N) + χ_ mom^ LO(N,)
    
    where the function χ_s(/M) is defined as the dual to the LO anomalous dimension γ_0,
    γ_0χ_s/M = M     ⇔    χ_s1/γ_0(N)=N,
    
    the function
    χ̃_0(M,N) = C_A/π[ψ(1) + ψ(1+N) - ψ(1+M) - ψ(2-M+N)]
    
    is the off-shell extension of the LO BFKL kernel with double counting with χ_s subtracted, and
    the function
    χ_ mom(N,) = c_ mom() f_ mom(N),     
    f_ mom(N) = 4N/(N+1)^2
    restores momentum conservation, i.e. the constraint γ_ DL(N=1,)=0 which translates into χ_Σ(0,1,)=1,
    through a suitable coefficient c_ mom. Because of the definition Eq. (<ref>), χ_s(/M) in M=0 equals 1,
    so we have
    
    c_ mom^ LO() = -χ_s/2 - χ̃_0(0,1).
    
    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. <cit.>, 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. <ref>.
    
    In order to obtain the coefficient of the -expansion of the DL anomalous dimension,
    we substitute the formal expansion
    γ_DL-LO(N,) = γ_0(N)+^2γ̃_1(N)+…
    into Eq. (<ref>) with χ_Σ given in Eq. (<ref>), and expand the equation in powers of .
    We stress that we are implicitly assuming that γ_0(N) in Eq. (<ref>) is the same used in the definition
    of χ_s, Eq. (<ref>); this has to be correct because of the way χ_Σ^ LO Eq. (<ref>) is constructed,
    as we shall verify shortly.
    On the other hand, γ̃_1 is a prediction of the resummation (hence the tilde).
    The tricky part for performing the expansion is the first (collinear) χ_s in Eq. (<ref>),
    as its argument /M is of (^0). For this term we then compute the expansion as
    
        χ_s/γ_DL-LO(N,)   = χ_s1/γ_0
            1-γ̃_1/γ_0+(^2)
           = χ_s1/γ_0 - γ̃_1/γ_0^2χ_s'1/γ_0 + (^2)
           = N + γ̃_1/γ_0' + (^2),
    
    where in the last equality we have used the definition Eq. (<ref>), assuming that γ_0 is the one appearing in Eq. (<ref>),
    and the formula for the derivative
    χ_s'1/γ_0 = -γ_0^2/γ_0',
    
    which can be obtained by deriving both sides of the first of Eq. (<ref>) with respect to /M.
    Here γ_0' denotes a derivative with respect to N.
    All the other terms can be simply expanded in powers of . Up to the first non-trivial order we get
    
    N = N +γ̃_1/γ_0' + χ_01/1+N + χ̃_0(0,N) - 
        χ_01/2 + χ̃_0(0,1)
    f_ mom(N) + (^2),
    
    where χ_01=χ_s'(0)=C_A/π, from which it immediately follows
    γ̃_1(N) = -γ_0'(N) 
        χ_01/1+N + χ̃_0(0,N)
         -χ_01/2 + χ̃_0(0,1)f_ mom(N)
    .
    
    Note that the (^0) term cancels automatically, which confirms that indeed the LO part of γ_DL-LO
    is given by the same γ_0 appearing in Eq. (<ref>).
    Now, it happens that, due to the explicit form of χ̃_0(M,N), Eq. (<ref>),
    χ̃_0(0,N) = C_A/π
        ψ(1+N)-ψ(2+N)
     = -C_A/π1/1+N,
    
    and hence we find
    γ̃_1(N) = 0.
    
    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. (<ref>)).
    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
    
        χ_Σ^ NLO(M,N,)
             = χ_s, NLO(M,) + χ_s, NLO(1-M+N,) 
           + χ̃_0(M,N)  + ^2 χ̃_1(M,N) + ^2χ_1^ corr(M,N) 
           + χ_ mom^ NLO(N,).
    
    Here, χ_s, NLO(M,) is the generalization of χ_s, constructed as the exact dual of the
    NLO anomalous dimension[As before, we use an approximate NLO anomalous dimension, see App. <ref>.]γ_0(N)+^2γ_1(N), which is an input at this order. This kernel satisfies the formal expansion
    χ_s, NLO(M,) = ∑_j=0^∞^j ∑_k=1^∞χ_jk/M^k;
    
    the j=0 term corresponds to χ_s(/M), Eq. (<ref>).
    The kernel χ̃_1(M,N) was given in Eqs. (A.23)–(A.29) of Ref. <cit.>, and we do not report it here.
    The extra term χ_1^ corr(M,N) takes into account running coupling corrections; its correct expression is
    χ_1^ corr(M,N) = β_0
        -C_A/πψ_1(2-M+N) - 1/(1-M+N)^2χ_s'/1-M+N + χ_0(M,N)/M-C_A/π M^2
    .
    
    This equation corrects Eq. (A.18) of Ref. <cit.> (i.e. Eq. (6.19) of Ref. <cit.>), 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, χ_ mom^ NLO(N,) restores momentum conservation, in the same form as Eq. (<ref>), with
    
    c_ mom^ NLO = -χ_s, NLO(2,) - χ̃_0(0,1) - ^2 χ̃_1(0,1) - ^2χ_1^ corr(0,1).
    
    Note that since χ_s, NLO(M,) is the exact dual of the NLO anomalous dimension, it equals 1 in M=0.
    Now we consider the expansion of the DL-NLO anomalous dimension
    γ_DL-NLO(N,) = γ_0(N)+^2γ_1(N)+^3γ̃_2(N) + …,
    
    where both γ_0(N) and γ_1(N) are assumed to be those used in the definition of χ_s, NLO
    (as before, this will be confirmed by the explicit computation),
    and γ̃_2(N) is what we aim to find.
    The expansion of χ_s, NLO is obtained by using the same technique used in Eq. (<ref>), and leads to
    χ_s, NLO(γ_DL-NLO(N,),) = N + ^2γ̃_2/γ_0' + (^3).
    
    Note that up to this order the expanded kernel χ_s(/M)+χ_ss(/M),
    corresponding to the first two terms in the j-sum of Eq. (<ref>) and originally used in ABF <cit.>,
    gives identical results.
    Substituting Eq. (<ref>) into Eq. (<ref>) with the NLO kernel Eq. (<ref>) and expanding
    in powers of  using Eq. (<ref>) we find the following expression
    
        γ̃_2(N) = -γ_0'(N)
             [ χ_11/1+N + χ_02/(1+N)^2 +χ̃_1(0,N) + χ_1^ corr(0,N) 
            -χ_11/2 + χ_02/4 + χ̃_1(0,1) +χ_1^ corr(0,1)f_ mom(N) 
           + C_A/π
            ψ_1(1+N)-ψ_1(1)
        γ_0(N) ].
    
    The coefficients of the expansion of χ_s, NLO are given by (see App. <ref>)
    
        χ_02   = -11C_A^2/12π^2 + n_f/6π^2(2C_F-C_A)
        χ_11   = -n_f/36π^223C_A-26C_F.
    
    Now, from the definition of χ̃_1 (see Ref. <cit.>)
    χ̃_1(M,N) = χ̃_1^ u(M,N) - χ̃_1^ u(0,N) + χ̃_1^ u(0,0),
    
    we immediately find
    χ̃_1(0,N) = χ̃_1^ u(0,0)
    
    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. <cit.>)
    χ̃_1(0,N) = 1/π^2
        -74/27C_A^2+11/6C_A^2ζ_2+5/2C_A^2ζ_3 + n_f4/27C_A+7/27C_F-1/3C_Fζ_2
    .
    
    On the other hand we have
    χ_1^ corr(0,N) = -β_0C_A/πζ_2,
    
    which is again N-independent.
    We can then rewrite Eq. (<ref>) as
    
        γ̃_2(N) = -γ_0'(N)
             [ ρ + χ_11/1+N + χ_02/(1+N)^2 -ρ + χ_11/2 + χ_02/4f_ mom(N) 
           + C_A/π
            ψ_1(1+N)-ψ_1(1)
        γ_0(N) ]
    
    with
    
        ρ   = 1/π^2
            -74/27C_A^2+11/12C_A^2ζ_2+5/2C_A^2ζ_3 + n_f4/27C_A+7/27C_F+1/6C_Aζ_2-1/3C_Fζ_2
         .
    
    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
    γ_0(N) = C_A/π N + (N^0),     γ_0'(N) = -C_A/π N^2 + (N^-1),     
    we have that up to NLL the singular behaviour of γ̃_2 in N=0 is given by (using f_ mom(0)=0)
    
        γ̃_2(N)    = C_A/π N^2
            ρ+χ_11+χ_02-2ζ_3C_A^2/π^2
        
           =C_A/π^3N^2C_A^2 (54 ζ_3 + 99 ζ_2 - 395) + n_f (C_A - 2 C_F) (18 ζ_2 - 71)/108,
    
    which indeed reproduces the correct NLL pole of the known three-loop anomalous dimension γ_2(N) <cit.>,
    while the LL 1/N^3 pole is accidentally zero.
    We stress that without the correction of the error in Eq. (<ref>) the constant in Eq. (<ref>)
    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. <ref> we have presented a new running coupling resummation and provided its  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 <cit.>, 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,) to the DL results violates momentum,
    which is further violated in the LO+LL^' by the matching function
    and in the NLO+NLL by Δγ_ss^ rc and its matching function.
    Momentum conservation can be restored by simply adding a function proportional to f_ mom(N), Eq. (<ref>).
    In summary, we have[Note that we are here using a notation for the resummed results, “res (N)LL^(')”,
    which differs from the one used in Ref. <cit.>, “(N)LO+(N)LL^(')”.
    The reason is that we will now use the latter name for the actual resummed results matched to any fixed order,
    while in these resummed results the fixed order is only approximate.]
        γ_+^ res LL (N,)
           = γ_DL-LO (N,) + Δ_DL-LOγ_ rc^  LL(N,) - Δ_DL-LOγ_ rc^  LL(1,) f_ mom(N), 
        γ_+^ res LL'(N,)
           = γ_DL-LO (N,) + Δ_DL-LOγ_ rc^ NLL(N,) + γ_ match^ LO+LL' (N,) 
              - 
            Δ_DL-LOγ_ rc^ NLL(1,) + γ_ match^ LO+LL' (1,)
         f_ mom(N), 
        γ_+^ res NLL(N,)
           = γ_DL-NLO(N,) + Δ_DL-NLOγ_ rc^ NLL(N,) + Δγ_ss^ rc(N,) + γ^ss_ match(N,) 
              - 
            Δ_DL-LOγ_ rc^ NLL(1,) + Δγ_ss^ rc(1,) + γ^ss_ match(1,)
         f_ mom(N),
    
    where the various functions have been introduced in Sect. <ref> and Sect. <ref>.
    Using the results in there, these expressions admit the following  expansions
    
        γ_+^ res LL(N,)     = γ_0(N) + ^2 β_0 3/32κ_0-c_01/N - f_ mom(N) + (^3),
        γ_+^ res LL'(N,)    = γ_0(N) + ^2 β_0 3/32κ_0-c_01/N - f_ mom(N) + (^3),
        γ_+^ res NLL(N,)     = γ_0(N) + ^2 γ_1(N) 
             + ^3 {β_0^2 κ_0/161/N - f_ mom(N)
                    -γ_0'(N)
          [ ρ + χ_11/1+N + χ_02/(1+N)^2 -ρ + χ_11/2 + χ_02/4f_ mom(N) 
                             + C_A/π
            ψ_1(1+N)-ψ_1(1)
        γ_0(N) ] }
        + (^4),
    
    with coefficients defined in Eqs. (<ref>), (<ref>) and (<ref>).
    
    Note that the functions γ_0 and γ_1 are those entering the definition of χ_s and χ_s, NLO
    in the DL kernel, which are not the exact fixed-order results (see App. <ref>).
    Thus, it is convenient to introduce the pure resummed contributions, defined by the difference
    between the above expressions and their expansion up to a given order, e.g.Δ_nγ_+^ NLL(N,) = γ_+^ res NLL(N,) - ∑_k=1^n ^k 
        γ_+^ res NLL(N,)
    _(^k),
    
    and similarly for the LL result. In this way, the resummed anomalous dimension matched to the (exact) fixed order is given by
    γ_+^N^nLO+N^kLL(N,) = γ_+^N^nLO(N,) + Δ_nγ_+^N^kLL(N,),
    
    where γ_+^N^nLO(N,) is the exact N^nLO anomalous dimension.
    In Ref. <cit.> we only considered the “natural” contributions
    
        Δ_1γ_+^ LL^(')(N,)     = γ_+^ res LL^(')(N,) - γ_0(N), 
        Δ_2γ_+^ NLL(N,)     = γ_+^ res NLL(N,) - γ_0(N) - ^2 γ_1(N)  .
    
    With the results of Eqs. (<ref>)
    we can now also compute
    
        Δ_2γ_+^ LL^(')(N,)     = γ_+^ res LL^(')(N,) - γ_0(N) - ^2 β_0 3/32κ_0-c_01/N - f_ mom(N), 
        Δ_3γ_+^ NLL(N,)
           = γ_+^ res NLL(N,) - γ_0(N) - ^2 γ_1(N)  
             -^3 {β_0^2 κ_0/161/N - f_ mom(N)
                    -γ_0'(N)
          [ ρ + χ_11/1+N + χ_02/(1+N)^2 -ρ + χ_11/2 + χ_02/4f_ mom(N) 
                             + C_A/π
            ψ_1(1+N)-ψ_1(1)
        γ_0(N) ] },
    
    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. <cit.> for all the details on its resummation, and we report here its
    expansion in powers of ,
    γ_qg^ NLL(N,) =  h_0 + ^2 h_1 γ_0^ LL'(N)
    + ^3 
        h_2 γ_0^ LL'(N)γ_0^ LL'(N) -β_0 + h_1 γ_1^ LL'(N)
     + (^4),
    
    where h_0=n_f/3π, h_1=n_f/3π5/3 and h_2=n_f/3π14/9 are numerical coefficients
    (already introduced in Sect. <ref>),
    γ_1^ LL'(N) is the (^2) term of γ^ LO+LL'(N,), which can be read off Eq. (<ref>),
    and γ_0^ LL'(N) is given (up to a factor ) by Eq. (B.10) of Ref. <cit.> which we report here for convenience
    (with the notations of App. <ref>)
    γ_0^ LL'(N) = a_11/N + a_10/N+1.
    
    Note that in Ref. <cit.> two forms for γ_qg^ NLL, differing by subleading terms, were considered
    and used to estimate an uncertainty of the resummation. Up to the order of Eq. (<ref>),
    both expressions give identical results, the difference starting at (^4).
    Using Eq. (<ref>) it is possible to construct the resummed contributions
    
        Δ_2γ_qg^ NLL(N,)    = γ_qg^ NLL(N,)-  h_0 - ^2 h_1 γ_0^ LL'(N), 
        Δ_3γ_qg^ NLL(N,)    = γ_qg^ NLL(N,)-  h_0 - ^2 h_1 γ_0^ LL'(N) 
             - ^3 
            h_2 γ_0^ LL'(N)γ_0^ LL'(N) -β_0 + h_1 γ_1^ LL'(N)
        ,
    
    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γ_qg^ NLL, as described in Ref. <cit.>.
    
    
    
    
    
      §.§.§ 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
    
    P_ij^N^nLO+N^kLL(x,) = P_ij^N^nLO(x,) + Δ_n P_ij^N^kLL(x,)
    
    with
    
        Δ_n P_gg^ NLL(x,)    = (1-x)^2 1-√(x)^4 
            Δ_n P_+^ NLL(x,) - C_F/C_AΔ_n P_qg^ NLL,nodamp(x,) - D
        
        Δ_n P_qg^ NLL(x,)    = (1-x)^2 1-√(x)^4 Δ_n P_qg^ NLL,nodamp(x,) 
        Δ_n P_gq^ NLL(x,)    = C_F/C_AΔ_n P_gg^ NLL(x,) 
        Δ_n P_qq^ NLL(x,)    = C_F/C_AΔ_n P_qg^ NLL(x,)
    
    (and similarly at LL)
    where Δ_n P_+^ NLL and Δ_n P_qg^ NLL,nodamp are the inverse Mellin of Δ_nγ_+^ NLL and Δ_nγ_qg^ NLL, respectively.
    The constant D is given by
    
    D = 1/d(1)∫_0^1 dx  x (1-x)^2 (1-√(x))^4 
        Δ_n P_+^ NLL(x,) +1- C_F/C_AΔ_n P_qg^ NLL,nodamp(x,)
    ,
    
    where d(N) is the Mellin transform of (1-x)^2 (1-√(x))^4.
    Note that, with respect to Ref. <cit.>, 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  and n_f.
      This ensures that these contributions start at (^3) and (^4) 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 , version  onwards.
      This code essentially computes the resummation of coefficient functions (and equally of the qg anomalous dimension) as described in Ref. <cit.>,
      and partly summarized in Sect. <ref>. The splitting function matrix is then constructed and momentum conservation imposed.
      The objects which are computed are Δ_2P_ig^ LL, Δ_3P_ig^ NLL (i=g,q),
      Δ_2K_hg, 2_a,g, 2_a,g(m) (a=2,L) and 2^ CC_a,g(m) (a=2,L,3)[Massive DIS coefficient functions are further sampled for various values of the quark massess.]
      on a n_f, , 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 , version  onwards, where
      the ,x grid is interpolated (cubicly), the quark components of the splitting, matching and coefficient functions are computed by colour-charge relations,
      and Δ_1P_ij^ LL, Δ_2P_ij^ NLL, Δ_1K_hg, 1_a,i and 1^ (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. (<ref>), (<ref>) and (<ref>), as well as the analogous for the coefficient functions.
      The latter was already done in the massless case in Ref. <cit.>; explicit results for the massive case are presented in App. <ref>.
      Explicit results for the splitting functions are presented in App. <ref>.
    
    We underline that the second step is the slowest, as the  code needs to compute several integrals for the varius grid points and the various
    functions to be resummed.
    The last step, performed by the  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  code, making the inclusion of small-x resummation very handy.
    As far as PDF evolution and construction of DIS observables is concerned, the code  has been included in the  code <cit.>,
    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 , we present here some representative results for the small-x resummation of splitting functions and DIS coefficient functions obtained with the techniques described in Ref. <cit.> 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 <cit.> 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. <ref>.
    
    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. <ref>.
    
    In Fig. <ref> 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  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
    (the fixed-order ones are the same in both  and Q_0 at these orders).
    The number of active flavours is n_f=4, and the value of the strong coupling is =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. <ref> can be compared directly to the ones presented in our previous paper <cit.>.
    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 <cit.>, as detailed in Sect. <ref>,
    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/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. <cit.>),
    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. <cit.>,
    namely symmetrizing the difference between our default construction of the γ_qg resummation which uses the evolution function in Eq. (<ref>),
    and what we obtain by switching to a simpler and formally equally valid version, Eq. (<ref>).
    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. (<ref>), 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. (<ref>) for the resummation of running coupling effects (our default)
    and the variant obtained using the simpler, yet equally valid Eq. (<ref>).
    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 (^3), 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 x^-1log x.
    These accidental zeros are not present beyond NNLO and so we expect the yet-unknown N^3LO splitting functions
    to significantly deteriorate the stability of the perturbative expansion because of their ^4 x^-1log^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^3LO.
    
    
    
    
    
    
     §.§ 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 <cit.>, but the input
    LL^' anomalous dimension used for computing the resummed coefficients did, as explained in Sect. <ref>.[Additionally, we changed the overall large-x damping, which is now uniformly chosen to be (1-x)^2(1-√(x))^4,
    as for the splitting functions, Eq. (<ref>).]
    The updated results are shown in Fig. <ref> 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 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. <ref> the analogous of Fig. <ref> 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 _a,i, with a=2,L and i=g,q, of which the functions _a,i
    defined in Sect. <ref> are the resummed contributions.
    For charm production (upper plots) we use =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 =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 n_f=4 scheme, Fig. <ref>, which instead assumed only coupling to light quarks,
    to obtain a complete prediction (see e.g. the decomposition Eq. (<ref>) 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. <ref> between fixed-order and resummed contributions
    is very similar to that of the massless results in Fig. <ref>.
    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 , 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. <ref>),
    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. <ref> we  plot the resummed contributions 2_2,g(m_c) and 2_L,g(m_c) in NC
    and 2^ CC_L,g(m_c), 2^ CC_2,g(m_c) and 2^ CC_3,g(m_c) in CC[For CC, we assume the production of a charm quark together with a massless anti-quark. This fixes the sign of the F_3 contribution.]
    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_a^c, a=2,L, which is entirely determined by the charm coefficient function.
    
    
    In the upper plots of Fig. <ref> 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. <ref> the ratio of the individual contributions to the resummed contribution
    to the structure function F_2, 2_2,g(m_c).
    We stress that if the axial contribution proportional to g_A^2 which remains after factoring out the g_V^2+g_A^2 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. <ref>, 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_A^2, 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 fixed-order 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^3LO, 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,  version .
    A fast interface to these results is available through a new version of its companion code,  version .
    Both codes are publicly available at the webpage https://www.ge.infn.it/ bonvini/hell/.
    The fast  code can be directly used to compute PDF evolution and DIS cross sections through the public code  <cit.>.
    
    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 <cit.> in the case of a PDF fit that includes DIS-only data.
    For the near future, we look forward to implementing other processes in , 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 , 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 <cit.>.
    
    
    We are always indebted to Richard Ball and Stefano Forte for many discussions and guidance on this topic. 
    We are particularly grateful to Tiziano Peraro and Luca Rottoli for discussions and collaboration in the early stages of this project.
    We also thank Giovanni Ridolfi and Juan Rojo for discussions and comments on our manuscript.
    The work of MB was partly supported by the Marie Skłodowska Curie grant HiPPiE@LHC.
    
    
    § OFF-SHELL DIS COEFFICIENT FUNCTIONS
    
    
    In this section we report the computation and the results for the off-shell coefficient functions needed for the small-x resummation of DIS observables.
    We focus on the contributions with an incoming gluon, which is off its mass-shell, as this is what enters in the resummation formula.
    The relevant diagrams are shown in Fig. <ref>.
    The off-shellness of the incoming gluon regulates the collinear divergence of the produced quark pair,
    so the case with massless quarks is just a finite limit of the case with massive quarks.
    Therefore, we start considering the most general case, in which the gluon converts to a
    massive quark (pair), with mass m_1, and then after interacting with the vector boson the quark changes mass, m_2
    (and of course the opposite setup, as in the right diagram of Fig. <ref>).
    This accounts for all possible type of interaction:
    
      * m_1=m_2 is the case in which the boson is either a photon or a Z (neutral current)
    
      * m_1≠ m_2 is the case in which the boson is a W (charged current).
    
    Additionally, we consider a generic coupling which includes vector and axial currents,
    even though we will see that the couplings will mostly factor out.
    
    Note that, in practice, the charged-current case is relevant only if one of the two quark is massive and the other massless.
    Indeed only charm, bottom and top masses cannot be neglected, but the top contribution to DIS is always negligible.
    Therefore, the only processes for which two different non-zero masses would be needed are c+W→ b and b+W→ c.
    But these processes are suppressed by the CKM matrix element V_cb∼ 4  10^-2,
    and by the phase-space restrictions due to the quark masses.
    
    Therefore, the only significant combinations in charged current will involve at most a single massive quark.
    These combinations are c+W→ s,d and s,d+W→ c,
    on top of the fully massless contributions u+W→ s,d and s,d+W→ u.
    
    
    
    
     §.§ Calculation of DIS off-shell partonic cross section
    
    
    In this section we summarize the computation of the DIS off-shell partonic cross section in the most general case
    in which the quarks in the final state have different masses and their coupling to the boson contains vector and axial components.
    We consider the process
    
    V^*q + g^*k→Q̅p_3+Q'p_4,
    
    where V^* is the generic off-shell vector boson, g^* is the off-shell gluon and Q and Q' are quarks.
    The final state quarks are on-shell, and their flavour is in general different (to cover the charged-current case), so their masses are p_3^2=m_1^2,  p_4^2=m_2^2.
    
    The two diagrams contributing to this process are depicted in Fig. <ref>.
    The invariant matrix element for the two diagrams is given by
    
    iℳ^μρ=-i g_s t^c u̅p_4[γ^μg_V+g_Aγ^5kγ^ρ- 2 p_3^ρ/t-m_1^2
      -γ^ρk-2p_4^ρ/u-m_2^2γ^μg_V+g_Aγ^5]vp_3,
    
    where g_s, g_V, g_A are the strong, vector and axial couplings respectively, and t=k-p_3^2=q-p_4^2 and u=k-p_4^2=q-p_3^2.
    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
    
        g_V^Z    =
          
            +1/2-2e_Qsin^2θ_ w   Q=u,c,t
        
            -1/2-2e_Qsin^2θ_ w   Q=d,s,b
          ,
           
          g_A^Z    =
          
            +1/2   Q=u,c,t
        
            -1/2   Q=d,s,b
          ,
    
    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
    
        g_V^W    = 1/√(2)V_QQ',    
          g_A^W    = -1/√(2)V_QQ',
    
    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_V^2→ 2g_V^Zg_V^γ and g_A^2→0
    (the combination g_Vg_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[A transverse vector p_⊥ is defined to have components only along directions orthogonal to the time and z directions,
    i.e. in components p_⊥=(0,p_x,p_y,0)=(0,𝐩,0), where we have also defined the two-vector 𝐩.
    Note that p_⊥^2=-𝐩^2.]q_⊥, k_⊥ and Δ_⊥ defined by
    
        q^μ   = z_1 p_1^μ+ q_⊥^μ, 
        
        k^μ   = z_2 p_2^μ- k_⊥^μ, 
        
        p_3^μ   =1-τ z_1 p_1^μ + τ̅ z_2 p_2^μ + q_⊥^μ-Δ_⊥^μ,
        
        p_4^μ   =τ z_1 p_1^μ + 1-τ̅ z_2 p_2^μ - k_⊥^μ-Δ_⊥^μ,
    
    where p_1=√(s)/21,0,0,1 and p_2=√(s)/21,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_V^2.
    
    We define the parton-level off-shell hadronic tensor as the squared of the matrix element Eq. (<ref>)
    averaged over the off-shell gluon polarizations <cit.> W_μν=-k_⊥^ρ k_⊥^σ/k_⊥^2ℳ_νσ^†ℳ_μρ.
    
    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[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.] W_μν=
    -g_μν W_1
    + k_μ k_ν W_2
    -iϵ_μνρσk^ρ q^σ W_3
    +q_μ q_ν W_4
    +k_μ q_ν +q_μ k_ν W_5
    +i k_μ q_ν -q_μ k_ν W_6.
    
    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
    
        A_2^2   =q_⊥^μ q_⊥^ν/Q^2 W_μν, 
        
        2A_2^2 -3A_L^2 = A_g^2   =-g^μν+q^μ q^ν/q^2 W_μν, 
        A_3^2   =-iϵ_μναβp_2^α q^β/4p_2 · q W^μν,
    
    where A_2^2, A_L^2 and A_3^2 are directly related to F_2, F_L and F_3, respectively.
    Using the definitions of the kinematic variables in the high-energy limit, Eqs. (<ref>),
    the squared amplitudes Eqs. (<ref>) can be rewritten as
    
        A_2^2   =8g_s^2[p_1 · kp_2 · q/p_1 · p_2]^2 
           ×{g_V^2+g_A^2/m_1^2-tm_2^2-u[1-1/q^2 k^21-2p_1· p_3p_2 · p_4/m_2^2-u-2p_1 · p_4p_2 · p_3/m_1^2-t^2]
                -g_A^2m_1+m_2^2/q^2-g_V^2m_1-m_2^2/q^21/m_1^2-tm_2^2-u},
        
        A_g^2   =8g_s^2 p_1· k^2/p_1 · p_2^2{g_V^2+g_A^2[p_2 · p_3^2+p_2 · p_4^2+m_1^2+m_2^2/2q^2p_2 · q^2/m_1^2-tm_2^2-u
                                   +m_1^2+m_2^2-2q^2/2k^2+m_1^2-m_2^2^2/2k^2 q^2p_2 · p_3/m_1^2-t-p_2 · p_4/m_2^2-u^2]
                -g_V^2-g_A^2[3 m_1 m_2/k^2p_2 · p_3/m_1^2-t-p_2 · p_4/m_2^2-u^2+m_1 m_2/q^2p_2 · q^2/m_1^2-tm_2^2-u]},
        
        A_3^2   = 4 g_s^2 2 g_V g_Ap_1 · k^2/p_1 · p_2^2
           ×{p_2 · p_3-p_2 · p_4/p_2 · q[q^2/k^2p_2· p_3/m_1^2-t-p_2 · p_4/m_2^2-u^2-p_2· q^3/m_1^2-tm_2^2-u]
                -m_1^2-m_2^2/k^2p_2 · p_3/m_1^2-t-p_2 · p_4/m_2^2-u^2}.
    
    These squared amplitudes must be integrated over the final state two-body phase space.
    In terms of the kinematic variables Eqs. (<ref>) we can express it as
    
        dΦ   =ν/2π^2 dτ dτ̅ d^2Δ δτ̅1-τν-𝐪-Δ^2-m_1^2 δτ1-τ̅ν-𝐤-Δ^2-m_2^2
           =1/8π^2dτ/τ1-τd^2Δ̃ δν-Δ̃^2+1-τm_2^2+τ m_1^2/τ1-τ-𝐪-𝐤^2,
    
    where we have further defined
    
        ν   =2z_1 z_2p_1· p_2,
        Δ̃   =Δ-τ𝐪-1-τ𝐤.
    
    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. (<ref>) by
    
        C_2(N,ξ,ξ_m_1,ξ_m_2)   =1/4π^2 g_V^2+g_A^2∫_0^1 dη η^N∫ dΦ A_2^2,
         C_L(N,ξ,ξ_m_1,ξ_m_2)   =1/4π^2 g_V^2+g_A^2∫_0^1 dη η^N∫ dΦ 1/3[2A_2^2- A_g^2],
         C_3(N,ξ,ξ_m_1,ξ_m_2)   =1/4π^2 2 g_V g_A∫_0^1 dη η^N∫ dΦ A_3^2,
    
    where we have introduced the variable
    η=Q^2/ν
    and we are expressing the result in terms of dimensionless ratios
    
        ξ   =𝐤^2/Q^2,     ξ_m=m^2/Q^2.
    
    Note that, as explained in the text, we are only interested in the N=0 Mellin moment.
    
    To carry out the integrations in Eqs. (<ref>) it is useful to express the following combinations
    
        m_1^2-t   =1/τ[1-τm_2^2+τ m_1^2+Δ̃-τ𝐤^2+τ1-τQ^2],
        
        m_2^2-u   =1/1-τ[1-τm_2^2+τ m_1^2+Δ̃+1-τ𝐤^2+τ1-τQ^2],
        τ̅   =[1-τ𝐪-𝐤-Δ̃]^2+m_1^2/1-τν,
        1-τ̅   =[Δ̃+τ𝐪-𝐤]^2+m_2^2/τν,
    
    in terms of the phase-space variables.
    In addition, it is convenient to use the following Feynman parametrizations
    
        1/m_1^2-tm_2^2-u   =∫_0^1 dyτ1-τ/[1-τm_2^2 +τ  m_1^2+Δ̃+y-τ𝐤^2+τ1-τ𝐪+y1-y𝐤^2]^2,
        1/m_1^2-t^2m_2^2-u^2   =∫_0^1 dy6τ^21-τ^2 y1-y/[1-τm_2^2 +τ m_1^2+Δ̃+y-τ𝐤^2+τ1-τ𝐪+y1-y𝐤^2]^4.
    
    In this way, most integrations are easy to perform,[Note that using the δ function in Eq. (<ref>) to perform the η integration imposes restrictions on the remaining integrations.
      In particular, the restriction can be cast as a lower integration limit for Δ̃^2.
      
      However, in the high-energy regime, this lower bound is immaterial, and the integral can be extended down to zero.]
    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
    
        x_1   =4y1-y,      x_2=4τ1-τ,
    
    and express the results as integrals over x_1 and x_2.
    General results for C_a0,ξ,ξ_m_1,ξ_m_2, 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.
    
    
    
    
    
     §.§ 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
    
        C_20,ξ,0,0   =/3π∫_0^1dx_1/√(1-x_1)∫_0^1dx_2/√(1-x_2)3/8ξ x_1+x_2^3
                   ×
            (2-x_1) x_2^2+ x_1 x_2 ξ (3 x_1+3 x_2-4 x_1 x_2)+(2-x_2)x_1^2 ξ^2
        
         C_L0,ξ,0,0   =/3π∫_0^1 dx_1/√(1-x_1)∫_0^1dx_2/√(1-x_2)x_2^2x_2+x_1ξ5-4x_1/4ξ x_1+x_2^3
         C_30,ξ,0,0   =0.
    
    These results coincide with those presented in Ref. <cit.>, 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. <cit.>.
    
    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
    
        C_20,ξ,ξ_m,ξ_m   =/3π∫_0^1dx_1/√(1-x_1)∫_0^1 dx_2/√(1-x_2)3/84ξ_m+ξ x_1+x_2^3
           ×{[2ξ^2 x_1^2-4ξ x_1^2x_2^2-ξ^2 x_1^2 x_2+3ξ x_1^2 x_2+3ξ x_1 x_2^2-x_1 x_2^2+2x_2^2
               +4ξ_m-ξ x_1^2 x_2+4ξ x_1-x_1 x_2^2-ξ x_1 x_2 -x_1 x_2+4x_2+16ξ_m^22-x_1 x_2]
               +8 g_A^2 ξ_m /g_V^2+g_A^24ξ_m+ξ x_1+ x_2^2},
        
         C_L0,ξ,ξ_m,ξ_m   =/3π∫_0^1dx_1/√(1-x_1)∫_0^1 dx_2/√(1-x_2)1/44ξ_m+ξ x_1+x_2^3
           ×{[x_2^2x_2+x_1ξ5-4 x_1+2ξ_m x_2ξ x_1^2-2ξ x_1-3 x_1 x_2 +4 x_2+8ξ_m^2 x_22-3 x_1]
               + 6 g_A^2 ξ_m/g_V^2+g_A^2[24ξ_m+ξ x_1^2+2+x_1 x_2^2+x_2ξ8-3x_1x_1+4ξ_m4+x_1]},
        
         C_30,ξ,ξ_m,ξ_m   =0.
    
    Note the presence in the above expressions of a term proportional to g_A^2/(g_V^2+g_A^2). 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 <cit.>,
    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. (<ref>) reduces to Eqs. (<ref>).
    
    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. (<ref>), 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
    
        C_20,ξ,ξ_m,0   =/3π∫_0^1 dτ∫_0^1 dx_1/√(1-x_1)3/24τξ_m+ξ x_1 +4τ1-τ^3
           ×{16ξ_m+1τ^2ξ_m+1-τ^2+ξ^2 x_1^2ξ_m+1-τ^2+τ^2
               +2ξ x_1^2 1-ττ3ξ_m^2+ξ_m6-16τ+16τ^2-16τ+3
               +8τ^2 x_1
            ξξ_m^2+ξ_m4-3τ+31-τ^2-ξ_m+1^21-τξ_m+1-τ
        },
        
         C_L0,ξ,ξ_m,0   =/3π∫_0^1 dτ∫_0^1dx_1/√(1-x_1)1/4τξ_m+ξ x_1 +4τ1-τ^3
           ×{16τ^2ξ_m+1-τ^23ξ_m+41-ττ+3ξ^2 x_1^2ξ_m
               +2ξ x_1^2 1-ττ9ξ_m^2+ξ_m9-32 τ-32τ1-τ
               +8τ^2 x_1
            ξ3ξ_m^2+10ξ_m1-τ+101-τ^2-3ξ_mξ_m+11-τξ_m+1-τ
        },
        
         C_30,ξ,ξ_m,0   =/3π∫_0^1 dτ∫_0^1 dx_1/√(1-x_1)3/24τξ_m+ξ x_1 +4τ1-τ^3
           ×{16τ^22τ-1ξ_m+1-τ^2+ξ x_1^2ξ2τ-1-61-ττξ_m+1-2τ
               +8τ^2 x_1
            1-τξ_m^2+ξ_m2-3τ+2τ^2-3τ+1+ξξ_m
        }.
    
    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_a0,ξ,0,ξ_m =  C_a0,ξ,ξ_m,0 for a=2,L,
    while there is a sign change in the parity-violating coefficient, C_30,ξ,0,ξ_m = - C_30,ξ,ξ_m,0
    (see also Ref. <cit.>).
    As expected, the massless limit of Eqs. (<ref>) reduces to Eqs. (<ref>).
    
    
    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:
    C̃_aN,M,ξ_m_1,ξ_m_2 = M ∫_0^∞ dξ ξ^M-1 C_aN,ξ,ξ_m_1,ξ_m_2,      a=2,L,3.
    
    In the massless case we find
    
        C̃_20,M,0,0   = /3π Γ^3(1-M)Γ^3(1+M)/Γ(2-2M)Γ(2+2M) 3(2+3 M-3 M^2)/2M(3-2M) ,
        C̃_L0,M,0,0   = /3π Γ^3(1-M)Γ^3(1+M)/Γ(2-2M)Γ(2+2M) 3(1-M)/3-2M , 
        C̃_30,M,0,0   = 0,
    
    which reproduce the results of Ref. <cit.>.
    In the massive case in neutral current we obtain
    
        C̃_20,M,ξ_m,ξ_m   =/3π 3/2ξ_m^MΓ^3(1-M ) Γ(1+M )/(3-2M)(1+2M)Γ(2-2M )
           ×{ 1+M-
            1+M-2+ 3 M -3 M^2/2ξ_m
         _2F_1 (1-M,1,3/2; -1/4ξ_m)
               -2g_A^2/g_A^2+g_V^21+2M/1-2M[8ξ_mM-2ξ_m-2
                                 +4ξ_m+1^2  _2F_11-M,1,-1/2,-1/4ξ_m]},
        
        C̃_L0,M,ξ_m,ξ_m   =/3π 3/2ξ_m^MΓ^3(1-M ) Γ(1+M )/(3-2M)(1+2M)Γ(2-2M )4ξ_m/1+4ξ_m
           ×{3+1-M/2ξ_m -
            3+1-M/ξ_m(1- M/4ξ_m)
         _2F_1 (1-M,1,3/2; -1/4ξ_m)
               -g_A^2/g_A^2+g_V^21+4ξ_m1+2M/12ξ_m^2[6ξ_m2M-3 _2F_12-M,1,3/2,-1/4ξ_m
                                           +3M-4 _2F_12-M,2,5/2,-1/4ξ_m]},
        
        C̃_30,M,ξ_m,ξ_m   =0.
    
    Apart from the contribution proportional to the axial coupling g_A^2, which is new,
    the other terms reproduce the results of Refs. <cit.>.
    Finally, the massive case in charged current is given by
    
        C̃_20,M,ξ_m,0   = /3π3/4ξ_m^M Γ1-M^3 Γ1+M/ M3-2M4M^2-1Γ2-2M1/ξ_m
           ×{MM^2ξ_m^2-2ξ_m-3-Mξ_m^2-4ξ_m-3+ξ_m+2
               +1-M1+ξ_mM^2ξ_m+3-3M-2 _2F_12M-1,1,M+1,1/1+ξ_m},
        
        C̃_L0,M,ξ_m,0   = /3π3/4ξ_m^M Γ1-M^3Γ1+M/M3-2M1+2M1-2MΓ2-2M1/ξ_m
           ×{2M^33ξ_m+1-M^2ξ_m^2+11 ξ_m+2-M ξ_mξ_m-2+2ξ_m^2 
                +M-12M^2ξ_m+1+Mξ_m^2-5ξ_m-2+2ξ_m^2
                       ×_2F_12M-1,1,M+1,1/1+ξ_m},
        
        C̃_30,M,ξ_m,0   =/3π3/2ξ_m^M1-MΓ1-M^3Γ1+M/M3-2MΓ2-2Mξ_m/1+ξ_m
           ×[_2F_12M,1,1+M,1/1+ξ_m-1-1/ξ_m].
    
    
     §.§ 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.
    
    
    
      §.§.§ Massless limit
    
    We can use the Mellin forms to study the massless limit ξ_m→0 of the resummed expressions.
    We have for neutral current
    
        lim_ξ_m → 0C̃_20,M,ξ_m,ξ_m   = C̃_2(0,M,0,0) + K̃_hg(M,ξ_m), 
        lim_ξ_m → 0C̃_L0,M,ξ_m,ξ_m   = C̃_L(0,M,0,0),
    
    and for charged current
    
        lim_ξ_m → 0C̃_20,M,ξ_m,0   = C̃_2(0,M,0,0)  + 1/2K̃_hg(M,ξ_m), 
        lim_ξ_m → 0C̃_L0,M,ξ_m,0   = C̃_L(0,M,0,0),
        lim_ξ_m → 0C̃_30,M,ξ_m,0   = C̃_3(0,M,0,0)  + 1/2K̃_hg(M,ξ_m), 
        lim_ξ_m → 0C̃_30,M,0,ξ_m   = C̃_3(0,M,0,0) - 1/2K̃_hg(M,ξ_m),
    
    where we have defined
    K̃_hg(M,ξ_m) = - /πξ_m^M1-M/MΓ^3(1-M)Γ(1+M)/(3-2M)Γ(2-2M).
    
    The function K̃_hg(M,ξ_m) contains the collinear singularity, appearing as a M=0 pole,
    and produces the logarithmic mass contributions when expanded in powers of M. In a sense, it represents the conversion
    from the collinear singularity regularized by the off-shellness and the very same singularity regularized by the mass,
    see Eq. (<ref>).
    
    
    The inverse Mellin of Eq. (<ref>), needed for the running coupling resummation Eq. (<ref>),
    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:
    
        Γ(1-M)Γ(1+M)    = M∫_0^∞ dt t^M-1/1+t, 
        
          4^1-M(1-M)Γ^2(1-M)/(3-2M)Γ(2-2M)   = ∫_0^1 dx x^1-M/√(1-x), 
        (4ξ_m)^M/M   = ∫_0^1 dy/y  (4yξ_m)^M.
    
    Then we change variable from t to ξ=4ξ_mty/x to write the product in the form of a Mellin transfrom:
    
        K̃_hg(M,ξ_m)
           = - /4π M∫_0^∞ dt t^M-1/1+t∫_0^1 dx x^1-M/√(1-x)∫_0^1 dy/y  (4yξ_m)^M 
           = - /4π M∫_0^∞ dξ ξ^M-1∫_0^1dx/√(1-x)∫_0^1dy/yxy4ξ_m/xξ+y4ξ_m.
    
    At this point we can simply read off the inverse Mellin transform in integral form
    
        𝒦_hg(ξ,ξ_m)
           = - /4π∫_0^1dx/√(1-x)∫_0^1dy/yxy4ξ_m/xξ+y4ξ_m
           = - /3π
            4ξ_m/ξ+logξ_m/ξ+ 2-4ξ_m/ξ√(1+4ξ_m/ξ)log√(ξ/4ξ_m)+√(1+ξ/4ξ_m)
        .
    
    which has been computed explicitly in the second line. The ξ-derivative of this expression is Eq. (<ref>).
    
    
    
    
      §.§.§ 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  expansion of the resummed results.
    For matching up to NNLO, i.e. (^2), we need the expansions up to (M).
    In the massless case we have
    
        C̃_20,M,0,0   =/3π[1/M+13/6+71/18-ζ_2M+M^2],
        C̃_L0,M,0,0   =/3π[1-M/3+M^2],
    
    while in the massive case in neutral current we obtain
    
        C̃_20,M,ξ_m,ξ_m   =/3π[1+4√(1+4ξ_m)^-12√(ξ_m)/2+g_A^2/g_A^2+g_V^212ξ_m 	^-12√(ξ_m)/√(1+4ξ_m)
                 +M/6√(1+4ξ_m){√(1+4ξ_m)3lnξ_m+5+2^-12√(ξ_m)
                   ×13-10ξ_m+61-ξ_mlnξ_m-61-ξ_m H_-,+-1/√(1+4ξ_m)
                   +12ξ_mg_A^2/g_A^2+g_V^2(8ln(4ξ_m)-16ln1+√(1+4ξ_m)
                     +^-12√(ξ_m)6lnξ_m+28-3H_-,+-1/√(1+4ξ_m))}
                 +M^2],
        
        C̃_L0,M,ξ_m,ξ_m   =/3π[√(1+4ξ_m)1+6ξ_m-8ξ_m1+3ξ_m^-12√(ξ_m)/1+4ξ_m√(1+4ξ_m)
                 +g_A^2/g_A^2+g_V^22ξ_m47ξ_m+2^-12√(ξ_m)-√(1+4ξ_m)/31+4ξ_m√(1+4ξ_m)
                 +M/31+4ξ_m√(1+4ξ_m){√(1+4ξ_m)12ξ_m-1+31+6ξ_mlnξ_m
                   +^-12√(ξ_m)6+81-6ξ_mξ_m-24ξ_m1+3ξ_mlnξ_m
                   +12ξ_m1+3ξ_m H_-,+-1/√(1+4ξ_m)
                   +2ξ_m g_A^2/g_V^2+g_A^2(2^-12√(ξ_m)23+82ξ_m+67ξ_m+2lnξ_m
                           -√(1+4ξ_m)8+3lnξ_m-6ξ_mln√(1+4ξ_m)-1/√(1+4ξ_m)+1
                           -62+7ξ_mH_-,+-1/√(1+4ξ_m))}
                 +M^2],
    
    having defined the harmonic polylogarithm
    
    H_-,+z=_21-z/2-_21+z/2+1/2ln1-z^2/4ln1-z/1+z.
    
    In the charged-current case, with only one massive quark, we have instead the following expansion:
    
        C̃_20,M,ξ_m,0   =/3π[1/2M+8-3lnξ_m+6ln1+ξ_m/6
                 +M/36(86-3lnξ_m10+3lnξ_m+6ln1+ξ_m13+3ln1+ξ_m
                       -36_21/1+ξ_m)+M^2],
        
        C̃_L0,M,ξ_m,0   =/3π[ξ_m/1+ξ_m1/2M
        +6+14ξ_m+3ξ_m1+2ξ_mlnξ_m-6ξ_m^2ln1+ξ_m/61+ξ_m
                 +M/361+ξ_m(-12+92ξ_m-18ξ_m^2ln^21+ξ_m+36ξ_m^2_21/1+ξ_m
                   +9ξ_m1+2ξ_mln^2ξ_m+6ξ_m7ξ_m-1lnξ_m
                   +66+ξ_m15-7ξ_mln1+ξ_m)+M^2],
        
        C̃_30,M,ξ_m,0   =/3π[-1/1+ξ_m1/2M-5+3+6ξ_mlnξ_m-6ξ_mln1+ξ_m/61+ξ_m
                 -M/361+ξ_m(56+301+2ξ_mlnξ_m+91+2ξ_mln^2ξ_m
                   -60ξ_mln1+ξ_m-18ξ_mln^21+ξ_m+36ξ_m_21/1+ξ_m)
                 +M^2].
    
    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. (<ref>),
    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 (M^0) terms in these expansions, after subtraction of massless collinear singularities as in Eqs. (<ref>) and (<ref>),
    reproduce the known () contributions <cit.>.
    
    Higher-order corrections exist both for neutral and charged currents (see e.g. <cit.>), however not in a form we could compare our expansion to.
    
    
    
    § DETAILS ON NUMERICAL IMPLEMENTATION
    
    
     §.§ The evolution function U
    
    
    A key ingredient of the formalism for resumming coefficient functions is the evolution function U(N,ξ), defined in Eq. (<ref>).
    As discussed in Ref. <cit.>, computing it exactly with the resummed anomalous dimension γ_+
    (specifically, with the LL^' anomalous dimension) requires integrating γ_+ over all values of  from 0 to ∞,
    which is numerically inconvenient.
    Therefore, following Refs. <cit.>, we use an approximate expression for the evolution factor,
    where the anomalous dimension is assumed to depend on  only at LO, with 1-loop running.
    This leads to the ABF evolution factor <cit.>
    U_ ABF(N,ξ) = (1+r(N,)logξ)^γ_+(N,)/r(N,)
    with
    
    r(N,) = -Q^2 d/dQ^2γ_+(N,)/γ_+(N,)
    = ^2 β_0 d/dγ_+(N,)/γ_+(N,).
    
    Note that the ratio r(N,) is such that the approximation reproduces the correct derivative of γ_+ in .
    However, this effect is strictly speaking beyond the formal accuracy we work with, so one could ignore
    it and replace
    
    r(N,)→β_0.
    
    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  is ignored
    and the evolution factor becomes simply
    
    U_ f.c.(N,) = ξ^γ_+(N,).
    
    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,  is computed at ξ Q^2 in U, Eq. (<ref>), 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
    ξ_0=exp-1/β_0,
    
    and ξ-integration must be limited to the region ξ≥ξ_0.
    In the fixed-coupling limit,  is frozen at its value in Q^2, so all values of ξ are in principle accessible, i.e. ξ_0f.c.=0.
    
    In Ref. <cit.> we derived the resummation formula which had originally the form (schematically)
    
    C_ res(N) = ∫_ξ_0^∞ dξ  C(0,ξ)  d/dξU(N,ξ).
    
    Then, for convenience in the numerical implementation, we integrated by parts to get
    
    C_ res(N) = -∫_ξ_0^∞ dξ d/dξ C(0,ξ)  U(N,ξ),
    
    where the boundary term at infinity vanishes thanks to C, and the boundary term in ξ_0
    is assumed to vanish because
    
    U(N,ξ_0) = 0.
    
    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. (<ref>) 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= c_0, with c_0 given in Eq. (<ref>).
    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. (<ref>).
    Indeed the anomalous dimension γ_+ that we use has a more complex structure than γ_s. Additionally,
    we use the approximation Eq. (<ref>), which is typically finite in ξ=ξ_0:
    
    U_ ABF(N,ξ_0) = (1-r(N,)/β_0)^γ_+(N,)/r(N,).
    
    Note that using the variant Eq. (<ref>)U_ ABF(N,ξ_0) is either 0 or ∞ depending
    on the sign of the real part of γ_+.
    This implies that the two formulations Eqs. (<ref>) and (<ref>) 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
    
    U(N,ξ) = U_ ABF(N,ξ)  D_higher-twist(ξ),
    
    with
    
    D_higher-twist(ξ) =
    
        1-logξ/logξ_0^1+1/β_0
           ξ<1 
    
    1    ξ>1.
    
    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 ,
    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. (<ref>) are indistinguishable from
    those obtained with the undamped function, which confirms that the results of Ref. <cit.>
    are unaffected (from the point of view of U).
    On the other hand, results obtained using Eq. (<ref>) are now identical to those obtained using Eq. (<ref>),
    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. (<ref>) or (<ref>). Therefore, we can first compute the inverse Mellin transform of
    U(N,ξ) Eq. (<ref>) as function of ξ, U(x,ξ), which we tabulate for various values of , x and ξ,
    and we subsequently use it to compute the ξ integration for each observable.
    In  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.
    
    
    
    
    
    
    
     §.§ 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.[In the charged-current case, when the mass of the quark before and after hitting the W is different,
    the form of X generalizes to X=1/(1+(m_1+m_2)^2/Q^2).]
    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. (<ref>).
    
    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. (<ref>), (<ref>) and (<ref>)
    by a term of the form
    Θ_N(N,X) = X^N/(1-X)N+1 = 1+(N)
    
    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
    Θ_x(x,X) = θ(X-x) 
        1-x/X^1/1-X
    .
    
    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. (<ref>).
    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. <ref>, speeding up the computation of resummed massive coefficient functions.
    This would not be possible using the N-space formulation, Eq. (<ref>), 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.
    
    
    
    
    
    
    
    
    
     §.§ 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 () and (^2) 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. (<ref>), (<ref>) can be added to the exact
    fixed-order anomalous dimension, restoring the correct fixed-order part.
    
    This approach was already considered in both Refs. <cit.>,
    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. <cit.>, which is a somewhat simplified
    version of the one originally proposed in Ref. <cit.>,
    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
    
        γ_0(N)    = a_11/N + a_10 +(N) 
        γ_1(N)    = a_22/N^2 + a_21/N + a_20 +(N)
    
    where a_22=0 (accidental zero), we propose the following approximate expression,
    γ(N) = a_1/N + a_0 - (a_1+a_0)2N/N+1,
    
    which is valid both at LO and NLO with appropriate coefficients. At NLO they are given by
    
        a_1    =  a_11 + ^2 a_21, 
        
          a_0    =  a_10 + ^2 a_20,
    
    and at LO one simply neglects the (^2) terms; the coefficients are given by
    
        a_11   = C_A/π, 
        
          a_21   = n_f26C_F-23C_A/36π^2, 
        
          a_10   = 
                   -11C_A + 2n_f(1-2C_FC_A+4C_F^2)/12π, 
        
          a_20   = 1/π^2
            1643/24 - 33/2ζ_2 - 18 ζ_3 + n_f4/9ζ_2-68/81 + n_f^213/2187
        .
    
    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 and a_10 at LO and both a_21 and a_20 at NLO
    provides an excellent approximation of the anomalous dimension in the small-N region.
    The last term in Eq. (<ref>) is a subleading (N) contribution
    which ensures momentum conservation γ(1)=0.
    At large N, Eq. (<ref>) behaves as a constant,
    γ(N) N→∞→ -2a_1-a_0.
    
    We verified that -2a_1-a_0>-1 for all values of ,n_f that we consider.
    In particular, the worst case is obtained at NLO with n_f=3, where the condition γ(N→∞)>-1 is satisfied for < 0.558.
    For all other values of n_f, and at LO, the  range of validity is larger.
    Clearly, this is more than enough, as for such values of  the perturbative hypothesis is lost.
    Indeed, in our current numerical implementation we never consider energy scales for which >0.35.
    
    Another very important advantage of the approximation Eq. (<ref>) is that its inverse function (i.e. the dual) can be computed analytically:
    
        χ_s(M)
           = a_1+a_0-M + √((M+a_1-a_0)^2 + 8a_1(a_1+a_0))/2(2a_1+a_0+M)
    
    (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. (<ref>) 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. (<ref>) for χ_s allows to have a single layer of zero-finding routines,
    improving significantly the numerical stability and the speed of the code.[A word of caution is needed in the choice of the branch of the square-root.
    Along the Mellin inversion integration path, which is the only place in N space where the DL anomalous dimension
    is computed, the standard branch of the square-root is suitable for the collinear χ_s. However,
    for the anti-collinear χ_s, a different branch is needed, where the cut is placed on the negative imaginary axis,
    which avoids crossing the cut during integration.]
    The expansion of χ_s in power of  is given by
    χ_s(M) = a_11/M + a_11a_10^2/M^2 + a_21/M + (^3),
    
    which implies, according to Eq. (<ref>),
    
        χ_01   = a_11 = C_A/π, 
        χ_02   = a_11a_10 = -11C_A^2 + 2n_f(C_A-2C_F)/12π^2, 
        χ_11   = a_21 = n_f26C_F-23C_A/36π^2.
    
    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.
    
    
    
    
    
     §.§ Inverse Mellin transforms
    
    
    Here we compute the various ingredients for the inverse Mellin transforms of Eqs. (<ref>)–(<ref>).
    Using the approximate form of γ_0 given in App. <ref>, we have
    
        γ_0 (N)    = a_11/N -(a_10+2a_11) +2(a_11+a_10)/N+1, 
        γ_0'(N)    = -a_11/N^2 -2(a_11+a_10)/(N+1)^2.
    
    We can also write the function Eq. (<ref>) as
    
    f_ mom(N) = 4/N+1 - 4/(N+1)^2.
    
    In Eq. (<ref>) a number of products appear. When a power of 1/N multiplies a power of 1/(N+1),
    it is always possible to write it as a sum of powers of individual poles in either N or N+1.
    Therefore, most terms can be computed by means of the following inverse Mellin transforms,
    
        M^-1
            1/N^k+1
           = (-1)^k/k!log^kx/x, 
         M^-1
            1/(N+1)^j+1
           = (-1)^j/j!log^jx,
    
    where the ^-1 symbol is a shorthand notation for representing the Mellin inversion,
    ^-1
        f(N)
     = ∫_1/2-i∞^1/2+i∞dN/2π i  x^-N f(N).
    
    Additionally, the function ψ_1(N+1), multiplied by powers of 1/N or 1/(N+1), appears.
    The computation of these inverse Mellin is trickier. We start from
     M^-1
        ψ_1(N+1)
     = log x/x-1.
    
    To compute inverse Mellin of ψ_1(N+1) with products of 1/N and 1/(N+1),
    we compute consecutive convolutions of Eq. (<ref>) with the inverse Mellin
    of a single power of 1/N and 1/(N+1), which are given respectively by 1/x and 1.
    Starting from
    
        M^-1
            1/Nψ_1(N+1)
           = ∫_x^1dz/xlog z/z-1 = _2(1-x)/x, 
         M^-1
            1/N+1ψ_1(N+1)
           = ∫_x^1dz/zlog z/z-1 = ζ_2 - _2(x) + 1/2log^2x-log(1-x) log x,
    
    we can easily compute successive integrals by just integrating these results (as functions of z) in dz/x or dz/z.
    The relevant results are
    
        M^-1
            1/N^2ψ_1(N+1)
           = 2[_3(x)-ζ_3]-[_2(x)+ζ_2]log x/x, 
         M^-1
            1/(N+1)^2ψ_1(N+1)
           = 2[_3(x)-ζ_3]-[_2(x)+ζ_2]log x-1/6log^3x, 
         M^-1
            1/N^3ψ_1(N+1)
           = 3[ζ_4-_4(x)]+[_3(x)+2ζ_3]log x + 1/2ζ_2log^2x/x, 
         M^-1
            1/(N+1)^3ψ_1(N+1)
           = 3[ζ_4-_4(x)]+[_3(x)+2ζ_3]log x + 1/2ζ_2log^2x + 1/24log^4x.
    
    Because the  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 (^3) contribution in Eq. (<ref>) 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
    
        M^-1
            1/N^2ψ_1(N+1)-ζ_2
           = -2ζ_3/x + 2 - log x +(x) 
           =  M^-1
            -2ζ_3/N + 2/N+1 + 1/(N+1)^2
         +(x)
         M^-1
            1/(N+1)^2ψ_1(N+1)-ζ_2
           = -2ζ_3 - 1/6log^3 x +(x) 
           =  M^-1
            -2ζ_3/N+1 + 1/(N+1)^4
         +(x)
         M^-1
            1/N^3ψ_1(N+1)-ζ_2
           = 2ζ_3/xlog x + 3ζ_4/x -3 + log x +(x) 
           =  M^-1
            -2ζ_3/N^2+3ζ_4/N - 3/N+1 - 1/(N+1)^2
         +(x)
         M^-1
            1/(N+1)^3ψ_1(N+1)-ζ_2
           = 3ζ_4 + 2ζ_3log x + 1/24log^4 x +(x) 
           =  M^-1
            3ζ_4/N+1-2ζ_3/(N+1)^2 + 1/(N+1)^5
         +(x)
    
    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. (<ref>).
    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.
    
    
    
    
    
    tocsectionReferencesjhep
