Small-$x$ phenomenology at the LHC and beyond: HELL 3.0 and the case of the Higgs cross section

Small-$x$ resummation has been proven recently to be a crucial ingredient for describing small-$x$ HERA data, and the inclusion of small-$x$ resummation in parton distribution function (PDF) determination has a sizeable effect on the PDFs even at the electroweak scale. In this work we explore the implications of small-$x$ resummation at the Large Hadron Collider (LHC) and at a Future Circular Collider (FCC). We construct the theoretical machinery for resumming physical inclusive observables at hadron colliders, and describe its implementation in the public code HELL 3.0. We focus on Higgs production in gluon fusion as a prototypical example, both because it is sensitive to small-$x$ gluons and because of its importance for the LHC physics programme. We find that adding small-$x$ resummation to the N$^3$LO Higgs production cross section can lead to an increase of up to 10% at FCC, while the effect is smaller (+1%) at LHC but still important to achieve a high level of precision.


Introduction
With the discovery of the Higgs boson, the Standard Model (SM) has been established as a successful theory of particle physics. While the SM cannot be the definitive theory, direct evidence of physics beyond the SM has not (yet) been observed at the LHC. The search for new phenomena beyond the SM at hadron colliders may be pursued by testing the SM to high precision, which is becoming possible thanks to the huge amount and excellent quality of the data collected by the LHC. To keep up, theoretical predictions must reach and possibly surpass the precision of the measurements. On the one hand, this requires refined theoretical predictions for the partonic cross sections for the processes of interest, which may be obtained by higher order computations, e.g. next-to-nextto-leading order (NNLO) or even next-to-next-to-next-to-leading order (N 3 LO) in some cases, and by the all-order resummation of important classes of logarithmic contributions. On the other hand, accurate and precise theoretical predictions for LHC processes require high-quality parton distribution functions (PDFs).
Recently, an important step forward towards improved determination of PDFs has been achieved in Refs. [1,2], where the resummation of small-x (high-energy) logarithms at next-to-leading logarithmic (NLL) accuracy as implemented in the public code HELL [3,4] has been included in PDF evolution and in the theoretical predictions of DIS observables. Small-x resummation has the important role of stabilizing the behaviour of DGLAP splitting functions at small x, which otherwise is compromised by powers of log 1 x . In particular, the first manifest instability appears at NNLO, and thus PDFs determined with NNLO theory are rather different to those determined with NNLO theory improved by NLL small-x resummation. This difference, determined at low Q 2 where the small-x HERA data lie, persists and is actually enlarged by DGLAP evolution at larger scales. As a result, resummed PDFs at the electroweak scale are very different from the NNLO ones at small x.
This raises an important question: how does this large effect impact LHC precision phenomenology? To properly answer, we need to compare fixed-order prediction with fixed-order (NNLO) PDFs to resummed predictions with resummed (NNLO+NLL) PDFs. While NNLO+NLL PDFs are now available, resummed predictions for LHC observables did not exist, or at least not in a format which makes them immediately usable for phenomenology. It is the goal of this paper to provide the theoretical setup to perform this resummation for inclusive observables with the public code HELL. The resummation of differential observables with HELL is left to future work.
As a first example of application of this setup, we will consider Higgs production in gluon fusion. Being initiated by two initial-state gluons, this process is very sensitive to the gluon PDF. Moreover, it is known that the inclusive Higgs cross section is dominated by contributions close to partonic threshold, which in turn implies that the gluon PDF contributes mostly at small x. In addition, the inclusive Higgs cross section in gluon fusion is known to N 3 LO [5][6][7][8], so we will provide all the ingredients to properly match small-x resummation of a physical process to N 3 LO for the first time. We then investigate the phenomenological implications of small-x resummation in Higgs production at the LHC, and to enlarge the sensitivity to the PDFs at small x also at higher-energy colliders, namely High-Energy LHC (HE-LHC) and a Future Circular hadron-hadron Collider (FCC-hh).
The structure of the paper is the following. In Sect. 2 we derive the formalism for smallx resummation of inclusive cross sections with two hadrons in the initial state. We discuss its implementation in the HELL code, and compare it to the original formulation [9] in the Altarelli-Ball-Forte (ABF) formalism [10][11][12][13][14][15]. We provide all the ingredients for matching small-x resummation in the partonic coefficient functions to N 3 LO. In Sect. 3 we move to Higgs production, and present first how the fixed-order cross section can be constructed to treat correctly the small-x behaviour at NNLO and N 3 LO, and then the effect of adding small-x resummation both at parton level and at the level of the physical cross section. We then draw our conclusions in Sect. 4, and collect technical details in App. A. This work represents a follow up of Refs. [3,4,16] and [1], and provides the foundations of Ref. [17].

Hadron-hadron collider processes at high-energy
The resummation of small-x logarithms in physical processes requires both using PDFs which include small-x resummation in their determination and evolution, and resumming to all orders the log 1 x contributions in the partonic coefficient functions. The latter resummation, which is the subject of this section, is based on the so-called k t factorization theorem, where the non-perturbative proton dynamics is factorized in parton distribution functions which depend on both the longitudinal momentum fraction x of the parton and its transverse momentum k t [18][19][20][21][22][23]. Relating this k tdependent PDFs to the usual collinear PDFs it is possible to resum the leading non-vanishing tower of small-x logarithms to all orders in the collinearly factorized partonic coefficient functions.
Another important ingredient for a stable small-x resummation is the inclusion to all orders of a class of subleading contributions originating from the running of the strong coupling α s [9,15]. In Ref. [3] the approach of Refs. [9,15] has been rederived and reformulated in a simpler and more general way, and proven to be identical to the original formulation under specific assumptions. The new formulation of Ref. [3] has been implemented in the computer code HELL [3,4], and it is very convenient from the analytical and numerical points of view, making the resummation of new processes and their inclusion in HELL rather straightforward. In Ref. [3], and subsequently in Ref. [4], this new formalism has been presented and used only for processes with a single hadron in the initial state, and specifically the deep inelastic scattering (DIS) process. In this section t H g g Figure 1. Leading order diagram for Higgs production in gluon fusion at hadron-hadron colliders. The quark running in the loop is predominantly a top.
we extend the formulation to processes with two hadrons in the initial state, relevant for hadronhadron colliders such as the LHC. This extension was already presented in the orignal formulation in Ref. [9,24]; in this section we will also show that our formulation, which is more general, reduces to the original one under the same assumptions considered for the single-hadron case.

Resummation formalism with two incoming gluon legs
We consider a hadron-collider process which is gluon-gluon initiated. The typical and cleanest example, which we will consider in greater detail later in Sect. 3, is Higgs production in gluon fusion, whose leading order diagram is depicted in Fig. 1. Other examples for which the results of this section will be relevant are, e.g., top-pair production and jet production.
We will assume that there are no collinear singularities to be subtracted at resummed level. Namely, the lowest order diagram with two gluons in the initial state must be finite without any collinear subtraction. Indeed, in order for a collinear singularity to be present, at least one of the gluons must split into a pair of quarks, one of which participates to the hard interaction. In other words, it must be possible to cut a quark line such that the diagram factorizes into a gluon splitting to quarks and a gluon-quark initiated subgraph. Therefore, in presence of collinear singularities in a gg initiated diagram, there must exist a lower order diagram which is gq initiated. But if this is the case, the resummation of the gg initiated process is subleading logarithmic with respect to the resummation of the gq initiated process, due to the extra power of α s and no logarithm in the g → qq splitting. Thus, at the leading non-vanishing logarithmic accuracy, contribution with two initial state gluons which require collinear subtractions do not contribute. This is the case for instance of Drell-Yan production, where indeed at lowest logarithmic order only the gq (and qq) channels contribute [24]. Of course, it is well possible that such gq channel contains itself a collinear singularity (as it happens in the Drell-Yan case). However, this process has a single gluon in the initial state, and the treatment is identical to the DIS case already discussed in Ref. [3].
Let us then focus on the cross section σ of a gluon-gluon initiated (at lowest order) process without collinear singularities, such as Higgs production, in hadron-hadron collision. In order to simplify the treatment, we take the Mellin transform of the cross section as where τ = Q 2 /s, with Q the hard scale of the process (e.g., the Higgs mass) and √ s the collider center-of-mass energy. The cross section in collinear factorization can be written in Mellin N space as the sum over partonic channels of simple products, where C ij are the collinearly factorized coefficient functions and f k the collinear PDFs, which depend on the factorization scale µ F ∼ Q. The strong coupling α s is in general evaluated at the renormalization scale µ R , which implies that there are logarithms of µ R /Q in the coefficient function to compensate its dependence; however, at the leading logarithmic accuracy we will consider, the µ R dependence is subleading, and we therefore omit it to simplify the notation. The factor σ 0 is chosen such that the coefficients functions are dimensionless, and normalized to 1 at LO in the dominant channel (supposed to be the gg channel in our case). In the high-energy limit, there is no need to distinguish the individual quarks, as they always contribute in the singlet combination. Thus, in this section, we will assume that the index q refers to the whole singlet PDF. 1 In the high-energy limit, the cross section can be also written according to the k t factorization theorem, which gives Here, F g is the k t -dependent gluon PDF, and C is the partonic coefficient function computed with two off-shell incoming gluons, the off-shellness being k 2 t . Obviously, the off-shell coefficient function is symmetric for the exchange of the two virtualities, k 2 t1 ↔ k 2 t2 . The k t -dependent gluon PDF can be related to the collinear PDFs through the relation [3,23,24] F g (N, where U is a function which is factorization scheme dependent. 2 In the Q 0 MS scheme [21,23,25,26] usually considered in the high-energy regime, and adopted also here, it is given by with [3] U (N, where γ + is the (small-x resummed) eigenvalue of the anomalous dimension singlet matrix which is singular at small x. With respect to Ref. [3], we are slightly changing the notation for the evolution function U , to extend it to the case in which µ F is different from the hard scale Q. More details on the actual form of the evolution function U and on the anomalous dimension used in its definition are given later in Sect. 2.2. Plugging Eq. (2.4) into Eq. (2.3) and comparing with Eq. (2.2), we find a relation between the coefficient functions in collinear factorization and the off-shell coefficient functions, With this assumption Eq. (2.2) is incomplete as it misses non-singlet contributions; however, this is irrelevant for the present discussion, which is focussed on the high-energy limit. 2 We observe that in Ref. [3] the same equation was written in terms of the "plus" eigenvector PDF, see Eq. (3.3) there. However, that expression misses a contribution +δ(k 2 t )f − (N, Q 2 ) in terms of the "minus" eigenvector PDF, which produces the δ(k 2 t ) term in Eq. (2.4). The results of Ref. [3] are unaffected; the only effect of that deficiency is that C − appearing in Eq. (3.26) could be written in terms of the off-shell coefficient function. However, that contribution is purely NLO, and could thus be extracted from the fixed-order computation.
where we have introduced the dimensionless variables Introducing the "auxiliary" coefficient function we can rewrite the quark coefficient functions as These expressions, already derived e.g. in Ref. [27], allow us to express both quark coefficient functions in terms of the gluon one and of the auxiliary function. Thus, from now on we will focus on the functions C gg , Eq. (2.7a), and C aux , Eq. (2.9).

The evolution function
The evolution function U , Eq. (2.6), is a key object for small-x resummation in partonic coefficient functions. Indeed Eqs. (2.7) encode small-x resummation thanks to the form of U , which contains the leading small-x logarithms to all orders, provided the anomalous dimension in there is itself accurate at least at LL. We thus now recall here some properties of this function already presented in Refs. [3,4], with particular focus on its µ F dependence that we are now including. First, we observe that the anomalous dimension in Eq. (2.6) is integrated between µ F and k t , and k t is integrated in Eqs. (2.7) over all accessible values. This means that the resummed anomalous dimension would be needed at all possible values of α s between zero and infinity, which represents a big numerical challenge. In order to avoid this problem, an approximation of the α s dependence of the anomalous dimension was proposed in Ref. [3], where (2.11) with β 0 the one-loop coefficient of the QCD β-function. Under this assumption, the evolution function becomes . (2.13) Note that α s in Eq. (2.13) is in principle α s (µ 2 F ); however, since the scale dependence of α s is subleading with respect to the leading logarithmic accuracy of the resummed coefficient functions, α s can be computed at any renormalization scale µ R without compensating for this change. The name ABF in Eq. (2.13) comes from the fact that with this approximated evolution function the approach of Refs. [9,15] is recovered, as proven in Ref. [3]. We will show in the next Sect. 2.3 that this is the case also for processes with two incoming gluons.
A second observation is related to the region of ξ accessible in the integrals Eqs. (2.7) and (2.9). As the strong coupling is running, the integration cannot extend beyond the position of the Landau pole Λ, identified by the equation where µ is in principle any scale. Solving the equation, we find that the smallest accessible value of ξ is where we have written ξ 0 both in terms of µ = Q and of µ = µ F . In particular, since the approximation Eq. (2.13) assumes α s to be computed at µ F , the last form is more adequate. Note that when ξ = ξ 0 the approximate evolution factor reduces to with α s = α s (µ 2 F ). This expression is in general finite; however, from general considerations (see Ref. [4]), we expect the evolution function to vanish in ξ 0 , at least at LL. The vanishing of U in ξ 0 is a property which turns out to be particularly useful, especially from a numerical point of view. Thus, to force the evolution function to vanish in ξ = ξ 0 , a non-perturbative higher-twist damping function was introduced in Ref. [4], (2.17) such that the final approximated expression for the evolution function is This expression, used throughout this paper, also allows to integrate by parts Eqs. (2.7) without producing any boundary term. Finally, we recall that in Ref. [3] a dedicated anomalous dimension, denoted LL , was constructed specifically for its usage in the evolution function U . This LL anomalous dimension is essentially a LL anomalous dimension, but its dominant small-x singularity is the one of the NLL result. However, in the recent work of Ref. [16] it has been suggested that this hybrid anomalous dimension may give rise to instabilities when expanded in powers of α s , as needed for the matching of resummed results to fixed order. Since the numerical limitations that lead to the introduction of the LL anomalous dimensions have been overcome in Ref. [4], it has thus been proposed in Ref. [16] to use directly the full NLL anomalous dimension, which also corresponds to the original approach of Ref. [15]. In this work we will consider both options later in Sect. 3, and we will provide further support to the suggestion of Ref. [16] of using the NLL anomalous dimension in the evolution function U . Thus, the new release of HELL, version 3.0, performs the resummation using the NLL anomalous dimension in U as default.
To conclude, we report the actual expressions that we will use for the resummation of coefficient functions, as implemented in the code HELL. On top of using the approximated evolution function Eq. (2.18), we integrate by parts so that the derivatives act on the off-shell coefficient function, and we compute the latter in N = 0, as its N dependence is subleading. The results are The second expression is equivalent to the result in the case of a single hadron in the initial state, such as DIS. The first equation is a new result. The actual numerical implementation of the first equation further requires (for numerical stability) a change of variables, as discussed in App. A.1.

Equivalence to the impact factor formulation
In this section, we show that Eq. (2.19a) leads formally to the same results as the formulation of Ref. [9]. The argument follows closely the one given in Sect. 3.3 of Ref. [3], extending it to the case of two initial gluons. We first introduce the so-called impact factor which is simply the double Mellin transform with respect to each k 2 t of the off-shell coefficient function. For later convenience, we have introduced in the definition of the impact factor a µ Fdependent term. Because by assumptions there are no collinear singularities, this function is analytic in M 1,2 = 0, and thus admits an expansioñ Because of the symmetry of the off-shell cross section for the exchange of the virtualities, the coefficients of this expansion are symmetric for the exchange of the indices,C kj =C jk . Now, we write again the off-shell cross section as the double inverse Mellin transform of Eq. (2.20), expanded as in Eq. (2.21), where we have used the identity We can now plug Eq. (2.22) into Eq. (2.19a) and get 3 This expression, computed at central scale µ F = Q, reproduces exactly the result of Ref. [9]. Indeed, the derivatives of the evolution function, due to its form Eq. (2.13), satisfy the recursion 4 which, together with the initial k = 0 condition U ht ABF (N, e ν ) ν=0 = 1, give rise to what is sometimes denoted γ k + with squared brakets [3,15,28]. In this notation the resummed result is written as which is a straightforward extension of the analogous resummation in the single-hadron case. Note that while using Eq. (2.26) is numerically challenging and necessarily approximate (the infinite series cannot be treated exactly in a numerical code), and its implementation cannot compete with the straightforward integral representation Eq. (2.19a), this form is quite useful for the expansion of the resummed result to fixed order, as we shall now see.

Expansion and matching to fixed order
The resummed results Eqs. (2.19), which contains the leading small-x contributions to all orders, are usually matched to a fixed-order contribution. To do so, we need to subtract from the resummed result its expansion in α s up to the fixed-order k considered, where the last sum is the truncated α s -expansion of the first (resummed) coefficient C. Then, this ∆ k C contribution is of O(α k+1 s ), and can be safely added to the fixed N k LO result. In this work, we consider the matching up to N 3 LO, which is the highest fixed-order accuracy available for Higgs production in gluon fusion. Thus, we need the expansion of the resummation up to O(α 3 s ). The construction of this expansion is obtained in a simple way using the impact-factor formulation, Eq. (2.26). To use it, we first write explicitly γ k + up to k = 3 (omitting the arguments for ease of notation), where r is given in Eq. (2.11). To proceed, we now need to expand in powers of α s both γ + and r. However, before doing so, we recall that in Refs. [3,4] a variant of the resummation, used to estimate the uncertainty from subleading contributions, was introduced in which r is replaced with α s β 0 , i.e. the α s -dependence of the anomalous dimension is treated as if it was just O(α s ), in line with the approximation Eq. (2.11). To cover both cases, up to N 3 LO it is sufficient to introduce a single parameter T , which equals 2 in the default case, and equals 1 in the limit r = α s β 0 . Introducing the expansion of the anomalous dimension With these expressions we can expand Eq. (2.28) as which can be now used in Eq. (2.26) to get the α s -expansion of the gg coefficient function: Depending on the anomalous dimension used in the evolution function U (see discussion in Sect. 2.2), which determines the actual form of γ 0,1,2 , this expression provides the first few orders of the resummed coefficient function needed to construct the resummed contribution ∆ k C gg up to k = 3.
To construct the expansion of the resummed coefficient functions for the other partonic channels, we need to expand the auxiliary function Eq. (2.19b). Straightforwardly, its impact-factor form can be derived from Eq. (2.26) by keeping only the j = 0 part of the sum, and flipping the sign Thus, its α s expansion is given by With these expressions it is then possible to construct also the resummed contributions ∆ k C qg and ∆ k C qq for the quark coefficient functions up to k = 3. All together, these expressions allow to match resummed results to N 3 LO. The computation of theC ij coefficients, needed for the expansions presented here, is detailed in App. A.2.

The first few orders of the anomalous dimension at LL and NLL
To conclude the section, we now present the analytic expressions of the O(α 1,2,3 s ) anomalous dimensions γ 0,1,2 needed for the matching of the resummed coefficient function to fixed order up to N 3 LO, Eqs. (2.32), (2.35) and (2.36). We treat both the case in which the anomalous dimension used is the LL introduced in Ref. [3] and the case in which the full NLL anomalous dimension is used, as suggested in Ref. [15,16], see discussion in Sect. 2.2. These expressions are obtained by expanding the purely resummed LL or NLL anomalous dimension, and have been already computed and presented in Refs. [3,4,16]. Thus, here we only report the final results [16]. For LL resummation we have while for NLL resummation the results are (2.42) The coefficients appearing above are and The two values of λ 1 , Eq. (2.44b), come from another variant of the resummation, used in the construction of γ + , which affects only the LL anomalous dimension at this order. More details can be found in App. A of Ref. [16]. All these expressions are implemented in HELL 3.0.

Resummed Higgs cross section at the LHC and beyond
We now turn our attention to a hadron-hadron collider process which is of great interest for LHC phenomenology: Higgs production in gluon fusion (ggH for short). Of course, Higgs physics is very interesting because the Higgs sector can be sensitive to new physics beyond the Standard Model. The inclusive Higgs cross section, which we are going to consider, is for instance sensitive to heavy particles coupling to gluons, which may then run in the loop of Fig. 1 and alter the production rate at the LHC.
Moreover, from a theoretical point of view, Higgs production is an interesting process because fixed-order perturbative QCD corrections are very large, with NLO being about twice the LO, and NNLO adding another ∼ 40% of the LO to the cross section. This motivated various studies to go beyond NNLO [28][29][30][31][32][33][34], culminating in the computation of this production process to N 3 LO [5][6][7][8] in the large top-mass limit. 5 It has been demonstrated in various ways that such large corrections mostly originate from soft-virtual contributions [28,30,33,34,38], dominant at large x, and can be resummed to all orders by means of threshold resummation techniques [29,[39][40][41], reaching N 3 LL accuracy [42][43][44][45][46]. 6 N 3 LO+N 3 LL predictions are very close to N 3 LO ones, thus suggesting that perturbative expansion is apparently converging, and giving some confidence that current theoretical predictions, such as the one recommended by the LHC Higgs Cross Section Working Group (HXSWG) [51], are sufficiently accurate for precision phenomenology.
In this work we investigate the effect of the all-order resummation of the small-x logarithms, i.e. those important in the opposite limit with respect to the one largely studied for this process. Indeed, Higgs production in gluon fusion is also one the LHC processes which is expected to be most sensitive to small-x logarithmic enhancement, due to the fact that it is gluon-gluon initiated at lowest order, and the gluon PDF is the most sensitive to small-x resummation effects. We will show that the consistent inclusion of small-x resummation has a sizeable effect. In the Q 0 MS scheme that we adopt, most of this effect comes from the use of resummed PDF instead of fixed-order ones, while the effect of resummation in the coefficient function is much milder. We conclude that including small-x resummation effects in Higgs production is essential to reach the accuracy advocated by the LHC HXSWG.
Going beyond LHC, the effect of small-x resummation is much larger at the High-Energy phase of the LHC, with √ s = 27 TeV, and even more at a Future Circular Collider (FCC) of √ s = 100 TeV. This is the case because increasing the collider energy, smaller and smaller values of x become accessible and increasingly important.
Before presenting resummed results, we recall that many results in the computation of the Higgs production cross section in gluon fusion are obtained within the so-called large top-mass effective field theory (EFT henceforth), where the top-quark is integrated out of the theory and its effect included as corrections in powers of m 2 H /m 2 t . However, in this theory the small-x region cannot be predicted correctly, as the x → 0 limit does not commute with the m t → ∞ limit of the EFT. Therefore, the correct inclusion of small-x resummation also requires a correct treatment of the small-x region at fixed order. In Sect. 3.1, we then first revisit how top mass dependence is included in fixed-order result and how the correct small-x logarithms can be included at fixed order. Then, in Sect. 3.2 and Sect. 3.3 we will show the impact of small-x resummation at parton level and on the cross section, respectively.
This section provides a detailed explanation of the small-x resummed results presented in Ref. [17] in the context of a double small-x plus large-x resummation.

Construction of the fixed-order result at small x with top mass dependence
The LO diagram for ggH production, Fig. 1, is a one-loop diagram with a massive internal particle. The NLO correction to this process has been carried out exactly [52,53]. However, from NNLO onward, the exact computation would require the evaluation of three-loop (or higher) diagrams with massive internal lines, which are out of reach of the current computational technology.
However, in the limit in which the partonic center-of-mass energy √ŝ is (much) smaller than twice the top mass m t , one can construct an effective field theory (EFT) in which the top loop shrinks to a point, leading to a pointlike interaction described by operators. The operator with the lower dimensionality does not depend explicitly on the top mass, except for a logarithmic dependence appearing in its Wilson coefficient. Operators with higher dimensionality will give rise to corrections suppressed by increasing powers of 1/m 2 t . Within this EFT it has been possible to push the computation of the ggH cross section at NNLO (both at the leading power level [54][55][56] and including a few corrections in 1/m 2 t [27,57,58]) and even at N 3 LO [5][6][7][8] (at leading power). Because the expansion parameter of the EFT iŝ where z = m 2 H /ŝ andŝ the partonic center-of-mass energy, the limits m t → ∞ and z → 0 do not commute. 7 Thus, the EFT cannot describe the small-z region correctly. For this reason, computations within the EFT can be (and have been) carried out as threshold expansions about z = 1, i.e. as power series in (1 − z), e.g. in Refs. [7,57]. Indeed, at large and medium z this expansion converges to the exact result, while at small z it is wrong anyway.
The goal of this subsection is to provide a way to supplement computations in the EFT with the exact small-z logarithms, which can be predicted from the resummation. Let's consider the generic coefficient function with perturbative expansion (omitting all unnecessary arguments and indices for simplicity) As we already stated, from NNLO onwards the exact m H /m t dependence is unknown. At NNLO, m H /m t effects have been computed as an expansion in up to the order p max = 3 [27,57,58], (3.4) while at N 3 LO only the first term is known (p max = 0) [5][6][7][8]. However, the expansion in m H /m t is not accurate at small z, since the actual expansion parameter is the one given in Eq. (3.1): in particular, the ρ t expansion is supposed to break down for z ρ t /4. Therefore, the small-z behaviour of the m H /m t expansion is unstable, exhibiting double-logarithmic enhancement and higher powers of 1/z at each extra order in ρ t , in contrast with the exact small-z behaviour which is single-logarithmic enhanced and always contains a single power of 1/z. The exact small-z behaviour, Eq. (3.6), can be predicted (at least at LL, i.e. n = k −1) from high-energy resummation. Our goal is therefore to understand how the exact Eq. (3.6) can be used to replace the wrong Eq. (3.5) of the large m t EFT computation. We recall that two different phenomenological solutions to this problem have been already proposed in Refs. [27,57] and [58], respectively. Here, we want to deal with this problem in a systematic way. As a first step, we need to understand whether the limit p max → ∞ converges to the exact result or not. At large z and for sufficiently small ρ t , the answer must be yes, or at least asymptotically yes. On the other hand, at small z the expansion clearly diverges, with new singularities appearing at each order in ρ t , Eq. (3.5). Thus, at small z, only the all-order sum may make sense, but certainly not any finite truncation of it. Therefore, in order to build up a sensible result, we need to make sure first to get rid of the bad small-z behaviour of the ρ t expansion, and once this is done we can add the exact small-z contribution, Eq. (3.6). The final expression must be such that the limit p max → ∞ tends to the exact result. 8 We will consider four possible approaches, in turn.

Method of subtraction
The first option that we consider consists in subtracting from the ρ t expansion the "wrong" small-z behaviour, Eq. (3.5), replacing it with the exact small z, Eq. (3.6). The resulting expression is where we have further introduced a function d(z), which represents a possible large-z damping to be uniformly applied to the small-z parts of the result. The role of this damping is to suppress the effect of the small-z contributions at large z: indeed, the 1/z terms without logarithms contained in the small-z contributions do not vanish at large z. This method, despite its simplicity and naturalness, has two important drawbacks. The first is that it requires the exact EFT result, and not just its (simpler to compute) threshold expansion, which, as we already mentioned, carries the same correct information, and differs only in the region where they are both wrong. At NNLO, the EFT small-z contribution is fully known for p = 0, 1, 2 [58], while the threshold expansion is also known for p = 3 [27,57]. At N 3 LO, only the leading term B (3) 0,0,5 was known until very recently, when in Ref. [8] the exact leading EFT result (p = 0) has been computed, thus allowing to use this method at N 3 LO as well. The second and perhaps more severe issue is that the result may strongly depend on the damping function used. Indeed, ideally, the two small-z contributions should cancel each other at large z. However, since the z → 1 limit of the small-z contribution in the ρ t expansion inherits its instability, there is no practical compensation at large z between what is subtracted and what is added. And this must not happen, since the C (k) p (z) terms are supposed to be reliable in the z → 1 limit. Thus the damping becomes a necessity, but its form is not prescribed by the procedure, leaving a degree of arbitrariness which may contaminate the result.

Method of threshold expansion
The expression in square brakets in Eq. (3.7) does no longer contain divergent terms at small z. 9 Thus, there is no loss of information if it is replaced with its threshold expansion, i.e. an expansion in powers of (1 − z). Let us define to be the function in square brackets in Eq. (3.7). Eq. (3.8) can be expanded at large z as 10 where a is a parameter, and the subscript "t.e." (threshold expansion) means that the function enclosed by those brackets is expanded in powers of 1 − z. In the equation above, the expansion coefficients c which is clearly not expandable in z = 1, and we have introduced the coefficients b j,n,i (a) according to which thus depend on the damping function d(z). The a parameter is in principle free, since the result is independent of a when the whole series in 1 − z is considered. Of course, any finite truncation of the series to i = i max will have a residual dependence on a, which can thus be used e.g. to estimate the uncertainty due to missing terms in the threshold expansion [7]. The coefficients C (k) p (z) have been computed in the first place as a threshold expansion at NNLO [27,57] and p,i are all known. Let us comment on the choice of the parameter a. The value a = −1 is the one adopted in Ref. [27,57] (as there the partonic cross section zC(z) is expanded). 11 This choice is such that terms behaving as 1/z are generated in the threshold expansion; however, these terms are not predicted correctly by the EFT, and indeed they have been subtracted in Eq. (3.8), so producing them can be dangerous. On the contrary, we note that if we choose a ≥ 0 both terms in Eq. (3.8) lead to a threshold expansion which does not grow at small z. In particular, for a = 0 the threshold expansion goes to a constant, while for larger a it vanishes in z = 0 (however a cannot be too large, otherwise it would affect the coefficient function in a region of medium z where the theshold expansion is reliable). This means that choosing a ≥ 0 the resulting coefficient function does not contain any leading small-z contribution. Thus, the threshold expansion with a ≥ 0 provides a natural and legitimate way of damping the EFT result at small-z.
This observation may suggest that, as long as the coefficient function is threshold-expanded with a ≥ 0, the term δC p (z), without subtracting the small-z terms. Indeed, at large z the two objects do not differ, due to the damping d(z) which suppresses the subtraction terms, and at small z both objects do not contain small-z contributions. Clearly, the subleading small-z terms (those not enhanced by 1/z) may differ, but these are anyway beyond our control, and certainly not predicted by the last term of Eq. (3.7). Thus, we may conclude that an equally good definition of the full result is 10 To simplify the following discussion, let us assume that for the gg channel the coefficient function is defined as the "regular" part of the decomposition C , where the distributional part contains plus distributions and δ(1 − z) functions, and the regular part everything else. 11 To be precise, in Ref. [27,57] also the distributional part in the gg channel is multiplied by 1/z, thus changing the actual c (k) p,i coefficients. However, this difference is immaterial for the present discussion.
provided a ≥ 0. In fact, Eq. (3.12) can be obtained with no approximations, by exploiting the dependence on the damping function d(z). Indeed, in this case, the damping function is no longer fully free, but we have a guiding principle how to choose its form. Namely, since at large z the first part of Eq. (3.12) is reliable up to O (1 − z) imax , the damping function must be suppressed at least as such that the exact small-z part does not spoil the accuracy of the threshold expansion. With this choice for d(z), the b j,n,i (a) coefficients are all zero for i ≤ i max : hence, up to i = i max , the EFT small-z terms do vanish. Thus, when using this damping function (or a more suppressed version of it), the threshold expansion of Eq. (3.7) gives exactly Eq. (3.12). Note that, because in Eq. (3.12) there is no subtraction of small-z EFT contributions, this formulation is simpler and very suitable for numerical implementation, both at NNLO and at N 3 LO.

Method of double subtraction
In Ref. [27,57] a different construction was considered, where the exact small z is added to the threshold expansion of C (k) p (z) after having subtracted from it its own threshold expansion, without applying any damping. To derive this possible approach, let us start again from Eq. (3.7), to which we replace the first part with its threshold expansion, and where we add and subtract the threshold expansion of the exact small z, (3.14) As far as the exact small-z term is concerned, it is clear that the damping function becomes unnecessary with this approach, as the term in square brackets in the last line of Eq. thus fixing the form of the b j,n,i (a) coefficients according to Eq. (3.11). The approach of Ref. [27,57] is obtained by ignoring the second and third terms of the first line of Eq. (3.14), such that the final result is We notice immediately that this form is very similar to Eq. (3.12), with the difference that the large-z behaviour of the exact small-z contribution is subtracted rather than being damped. The result is in both cases a small-z contribution which starts at O (1 − z) imax+1 , and with the same small-z behaviour: it should thus give similar results. This observation can be considered as an a posteriori argument to justify neglecting the second and third terms of Eq. (3.14). The sum of these terms isn't necessarily small, and in fact for finite p max it may be sizeable. The theoretical argument behind neglecting them could be that in the p max → ∞ limit they vanish. However, the limit is divergent, making such an argument hard to prove.

Method of generalized expansion
The method of threshold expansion, Eq. (3.12), can be straightforwardly generalized by expanding about a generic z = z 0 , log n z z , (3.17) where this time the damping function must be of O (z 0 − z) imax+1 , in order to avoid spoiling the accuracy of the expansion. Such a function can be where the functional form is such that in z → 0 the damping is ineffective, and the theta function ensures that the small-z contribution remains zero for values of z > z 0 where it is forced to vanish. Eqs. (3.17) and (3.18) reproduce exactly Eqs. (3.12) and (3.13) for z 0 = 1. 12 For z 0 < 1, the expansion about z 0 cannot be valid in a region close to z = 1, essentially because of the presence of logarithmic terms diverging in z = 1 which force the convergence radius to be strictly less than 1 − z 0 . However, this limitation can be simply overcome by patching this result with a purely threshold-expanded one at some z = z 1 , with z 0 ≤ z 1 < 1, to be used for all z > z 1 . The advantage of this approach is that the EFT result is used in an extended region of z, while the contribution from the exact small-z terms, which are only known at LL, is relegated to a region of smaller z. The limitation of this approach is that z 0 cannot be too small. Indeed, the EFT approach breaks down for z ρ t /4, so z 0 must be sufficiently larger than this value. For the physical Higgs and top masses, ρ t /4 0.13. An interesting value to consider is z 0 = 1/2, for two reasons. The first is that the EFT expansion parameter, ρ t /(4z), equals 0.26 in z = z 0 = 1/2, which is just twice as large as its value at threshold z = 1, and thus hopefully still sufficiently small for the EFT to be reliable. The second, more practical, is that the expansions coefficients of the leading EFT contribution (p = 0) have been computed for z 0 = 1/2 in the recent work of Ref. [8] up to N 3 LO, making the implementation of this method rather straightforward.

Conclusion
The considerations above bring us to the conclusion that the method of threshold expansion, Eq. (3.12), using the damping function Eq. (3.13) and a = 0 provides the best way of implementing the correct small-z logarithms in a EFT result, such as the NNLO and N 3 LO ones. We have implemented this method in the public code ggHiggs, version 4.0 onwards. The method of double subtraction, Eq. (3.16), has also been implemented in the code to test the sensitivity of the results on the method of including small-z contributions (this method was already used in previous versions of ggHiggs for the NNLO, with a = −1, as prescribed in Ref. [27,57]). The method of generalized expansion, Eq. (3.17), with z 0 = 1/2 has been implemented at N 3 LO to investigate the possibility of improving the description of the transition region 10 −2 z 10 −1 , as we will discuss below.
The actual numerical implementation of the exact small-z logarithms has to face with the limitation that we know from resummation only the leading contribution, A (k) n with n = k − 1, while the coefficients with n < k − 1 are unknown at NNLO and N 3 LO. Here we can follow two possible approaches: we can either include only the known A (k) k−1 , setting to zero all the other coefficients, or we can guess their values. Since the subleading logarithmic contributions are likely more important than the leading one in a region of medium-small z, keeping these coefficients certainly helps. However, exactly because they may be relevant, their values must be guessed wisely. 12 To be precise, for z 0 = 1, the coefficientsc To do so, we follow the idea proposed in Ref. [28], namely we include subleading contributions as predicted by the LL resummation. To be precise, we use Eqs. (2.32), (2.35) and (2.36) to construct the O(α 2 s ) and O(α 3 s ) contributions to the coefficient functions in the various partonic channels in terms of the coefficientsC ij . The anomalous dimensions γ 0,1,2 appearing in those equations are taken to be the expansion terms of the exact "plus" eigenvalue of the singlet anomalous dimension matrix, rather than the ones predicted by the resummation, which is appropriate for a fixedorder prediction. Finally, these expressions are expanded about N = 0 to identify the resulting A (k) n coefficients. This procedure is certainly not fully correct. However, in a NLL (or higher) resummation, these contributions will be part of the full prediction, which will also contain some additional correction due e.g. to the impact factor at the next perturbative order. The hope is that these corrections be less important that the contributions that we include, such that the prediction is at least reasonable.
It is clear, however, that such an implementation is not fully satisfactory, especially at N 3 LO, where only one out of three parameters is known, i.e. A are only guessed. Therefore, it is important to assess, at least qualitatively, the effect of not knowing all the small-x contributions. To do so, we can consider variations of the unknown parameters. There are various ways in which these could be done, none of them being particularly justified. Thus, we consider only a very simple variation, which is obtained by setting to zero the coefficient of 1/z, A (k) 0 . At NNLO, this is the only unknown coefficient, while at N 3 LO it is the one that governs the contribution which has the largest impact at medium z, and it is thus sufficient to infer an uncertainty for those values of z relevant for LHC or FCC. Incidentally, the resulting uncertainty covers the difference between the two implementations Eq. (3.12) and Eq. (3.16), which is located in the medium-z region, as well as the negligible effect of changing a from 0 to 1. Thus, we can take the uncertainty band obtained setting A (k) 0 = 0 as a good representative of all the small-z uncertainties at fixed order.
In Fig. 2 we show the partonic coefficient functions for the gluon-gluon, gluon-quark and quarkantiquark partonic channels for factorization and renormalization scales both equal to half the Higgs mass (m H = 125 GeV), which is nowadays the default central scale adopted by most groups [51], and with m t = 173 GeV. In the gg case, only the regular part of the coefficient function is shown, as defined in footnote 10. Results are presented in solid red (NLO in the upper plots, NNLO in the mid plots, N 3 LO in the lower plots). At NLO the result is exact [52], while beyond NLO it is constructed according to Eq. (3.12) with a = 0 and with damping function Eq. (3.13). Consequently, NNLO and N 3 LO results are supplemented with an uncertainty band, obtained as described above by setting A (2) 0 = 0 and A (3) 0 = 0 respectively, and symmetrizing the variation with respect to our central prediction. For each plot, the leading EFT approximation (p = 0) is also shown in dashed black. At N 3 LO, the alternative implementation Eq. (3.17) is also shown, together with its own uncertainty band, in dot-dashed blue. Note that at N 3 LO the small-x contributions are different depending on whether the MS or the Q 0 MS scheme is used. Here we decide to use the Q 0 MS scheme also at fixed order, to match the scheme adopted in the resummed results that we will consider in the next subsection.
Several comments are in order. First, it is apparent that the EFT approximation has the wrong small-z behaviour, as it exhibits double logarithmic enhancement (at leading power) rather than the correct single logarithmic enhancement. Indeed, as discussed before, the EFT is expected to fail for z ρ t /4 0.13: this is apparent from the plots, where the red and black curves behave differently for values of z smaller than about ρ t /4. An exception is the qq channel at NLO, where the contribution from the produced s-channel gluon is resonant at the tt threshold in z = ρ t /4, producing the peak which is clearly not present in the EFT approximation. In this case, the agreement between the EFT and exact result is restricted to a region of larger z. This effect is expected to be diluted at higher orders, due to the richer dynamics; nevertheless it also suggests 0.05, thus also reducing the impact of the uncertainty from subleading logarithms in the medium-z region. However, the considerations above suggest that z ∼ 0.05 is dangerously outside the region of reliability of the EFT (which is roughly speaking z > 0.2), so the gain in precision (smaller uncertainty) of this construction is compensated by a loss in accuracy (the unknown exact result may lie outside the estimated uncertainty). This suggests to discard the construction based on Eq. (3.17), and use the safer construction based on Eq. (3.12).
Finally, at NNLO we observe a reduction of the uncertainty band when going from gg to qg and to the purely quark initiated channel. This reflects a relatively less important contribution of the small-z logarithms in quark channels. At N 3 LO the pattern is the same, but the uncertainty bands are bigger, as appropriate due to the fact that the fraction of known small-z terms at this order is smaller. We stress that, in general, the displayed uncertainty is likely an overestimate of the actual uncertainty, since the coefficient A (k) 0 is brutally set to zero rather than varied in a reasonable range. Thus, the uncertainty band will be useful only to visualize the potential impact of subleading logarithmic contributions and to motivate further work towards their computation, rather than for computing an actual uncertainty on the cross section.

Impact of high-energy resummation at parton level
Having described how the exact small-z behaviour is included in fixed-order computations performed within the large top-mass EFT framework, we now investigate the effect of supplementing the fixedorder computation with the all-order resummation of small-z logarithms. At parton level, this is implemented by adding to the fixed-order coefficient functions the resummed contributions ∆ k C ij defined in Sect. 2.4. In this section we study the impact of resummation on partonic coefficient functions, while the effect on the physical cross section will be discussed in the next section.
In Fig. 3 we report the partonic coefficient functions in the same format as Fig. 2. Results are presented at fixed order in solid red, and with resummation (in the Q 0 MS scheme) in the two implementations: using the NLL anomalous dimension in dashed blue, and using the LL one in dot-dot-dashed yellow (see discussion in Sect. 2.2 and in Ref. [16]). The fixed-order results are also supplemented by the band which represents a rough estimate of the potential impact of unknown subleading logarithmic contributions, as described in the previous subsection. Similarly, the resummed results are supplemented by an uncertainty band, obtained varying subleading logarithmic contributions related to running coupling effects in the resummation procedure, as described in Refs. [4,16]. 13 At NLO and NNLO the two implementations of the resummation give qualitatively similar results, deviating from the fixed order for z 10 −1 at NLO and z 10 −2 ÷ 10 −3 at NNLO. The growth of the resummed contribution is slightly stronger when the NLL anomalous dimension is used. The uncertainty band of the LL variant is slightly larger than the one of the NLL variant, and covers the latter result everywhere, making them fully compatible. At NNLO, we notice that the resummed result lies within the fixed-order uncertainty band for z 10 −4 ÷ 10 −3 , which is the region most important for phenomenology. If the bands represented faithfully the uncertainty from unknown subleading logarithms at fixed order, then the effect of resummation in the partonic coefficient functions would be irrelevant compared to such uncertainty.
At N 3 LO the general pattern is similar, with some important differences. The resummed contribution, computed using the NLL anomalous dimension, is a small correction which lies entirely within the fixed-order uncertainty band for the whole z range shown. However, this time the behaviour of the resummed result with LL anomalous dimension is rather different. In general, the effect is larger than the corresponding one with NLL anomalous dimension, and no longer fully compatible with it, even though the uncertainty band is also increased. Moreover, in the gg and qg channels, there is a sizeable contribution of the resummation in a region of medium-large z, 10 −2 z 0.2. This is entirely due to the O(α 3 s ) expansion of the LL anomalous dimension γ LL 2 , Eq. (2.39), and is indeed absent in the quark-quark channel which does not depend on it, see Eq. (2.36). This large contribution is in a region of z which cannot be considered to be dominated by small-z logarithms, and therefore has to be interpreted as a spurious effect. Indeed, the all-order resummed results with NLL and LL anomalous dimensions agreed in that z region when matched to NLO and NNLO, so there is no physical underlying reason for which they should give such different results when matched to one order higher.  Our interpretation of the origin of this spurious behaviour is the fact that while the LL anomalous dimension makes perfect sense to all orders, its α s expansion may be unstable order by order, perhaps due to its hybrid nature. This is not the case for the NLL anomalous dimension, which has a well behaved α s expansion. This conclusion is in agreement with the analysis of Ref. [16], and represents another motivation for favouring the use of the NLL anomalous dimension in place of the LL one in the computation of resummed coefficient functions, in particular when these are matched to N 3 LO or to a higher order. Nevertheless, we must bear in mind that the resummed result based on the LL anomalous dimension differs from the NLL one by subleading contributions. Therefore, the difference between the two formulations probes unknown subleading logarithms. We see that this difference is similar (slightly more conservative) than the uncertainty band on the NLL-based resummed result when matched to NLO or NNLO, and could thus be used as an alternative way of estimating subleading logarithmic uncertainty. When matched to N 3 LO, this difference is rather larger than the simple blue band, especially in the medium-large z region, and using it as a subleading logarithmic uncertainty would be rather conservative. However, given that we do not really know how large these subleading logarithms may actually be, we suggest to use this difference as  a measure of such uncertainty. As we will see in the next section the resulting uncertainty at the physical cross section level is very reasonable.
Another powerful way of visualizing the effect of resummation at parton level is through the Mellin transform of the coefficient functions. In Fig. 4 we show the dominant one, C gg , as a function of the Mellin variable N for positive real N . This time we include the full coefficient function, and not just the regular part, since the Mellin transform of a distribution is an ordinary function. In fact, the distributional part of the coefficient function is responsible for the growth of the coefficient function at large N [59]. Moreover, it is known [38,59,60] that a saddle point evaluation of the Mellin inversion integral defining the full cross section (i.e. including both the coefficient functions and the PDFs) provides an excellent approximation to the exact result, thus showing that the bulk of the contribution of the coefficient function to the cross section is encoded in its value at the saddle point N = N saddle . From Ref. [38] we know that the saddle point for SM Higgs production varies from 14 N saddle 1.1 for LHC at √ s = 7 TeV to N saddle 0.9 for LHC at √ s = 14 TeV and to N saddle 0.7 for FCC at √ s = 100 TeV. Thus, the region of interest for phenomenology in a vast range of hadron-hadron collider energies is all located in a small range of N close to N = 1.
In the left plot of Fig. 4 the full LO, NLO, NNLO and N 3 LO coefficient functions are shown, together with the resummed NLO+LL, NNLO+LL and N 3 LO+LL counterparts. Since the plot becomes busy in the small-N region, we also plot in the right panel the ratio of each curve to the highest order curve, N 3 LO+LL. We see that the resummed results depart from the fixed order for N < 1, and they all diverge at the same N = N pole > 0, which is determined by the resummation. Thus, they all grow stronger than each fixed order, which instead are singular in N = 0. Interest-ingly, the N 3 LO+LL curve is very close to the N 3 LO curve even at rather small N 0.2, which is in line with the behaviour found in the z-space plots, and shows that the effect of small-z resummation on the N 3 LO coefficient function is expected to be negligible, since the saddle point is in a region where N 3 LO and N 3 LO+LL are almost identical. In particular, the effect of subleading logarithms at fixed order, estimated by the coloured filled bands in the plots, is likely more significant than the effect of all-order resummation, both at NNLO and N 3 LO. The fixed-order uncertainty bands also appear to be larger than the uncertainty on the resummed contributions, estimated as the difference between the NLL and LL variants, and shown with a pattern. While we may hope, as already discussed, that these bands be over conservative, it seems important to take this observation as a strong motivation to work towards improving the knowledge of the small-z behaviour of the Higgs partonic coefficient functions.

Impact of high-energy resummation on the cross section
We now move to the physical cross section. It is defined as the convolution of the partonic coefficient functions with the PDFs, according to Eq. (2.2) which in momentum space reads 15 (3.20) with τ = m 2 H /s and s che collider center-of-mass energy, and we have restored the dependence on the renormalization scale µ R . Since high-energy resummation affects PDF evolution, and PDFs at smallx are mostly determined by HERA data at low Q 2 which are thus very sensitive to resummation effects and very "far" from the Higgs scale, it is crucial to use PDFs which have been determined and evolved using resummed theory when computing physical predictions which include high-energy resummation.
Recently, such PDFs have been determined in the context of the NNPDF methodology to PDF fitting [1]. Soon after, the xFitter collaboration also performed an analogous determination [2], whose findings are in agreement with those of the NNPDF study. In both cases, PDF sets have been fitted using fixed-order theory (NLO or NNLO 16 ) supplemented by high-energy resummation at NLL in the Q 0 MS scheme provided by the HELL code, version 2.0. To be precise, resummation in DGLAP evolution is really NLL, while resummation in DIS coefficient functions is just formally NLL, since the LL contribution vanishes. In this case, we would refer to the accuracy of resummation in DIS as absolute NLL but relative LL (for this notation, see Ref. [4]). In this respect, Higgs resummation, which is relative LL, is consistent with the PDF sets of Refs. [1,2].
In fact, since the Higgs cross section is known at fixed order up to N 3 LO, a consistent computation would require the use of PDFs obtained with N 3 LO theory, supplemented by resummation when computing resummed cross sections. However, this would require four-loop DGLAP splitting functions, which are not known yet, even though recently there has been some impressive progress towards their computation [35][36][37]. Therefore, for the time being we can only rely on NNLO (or NNLO+NLL) PDFs.
We will focus on the PDFs of Ref. [1], which are publicly available. In that work, various families of PDF sets have been obtained by using different datasets. The mainstream family is based on a global dataset, which includes on top of DIS data a large amount of "hadronic" data (mostly Drell-Yan, jet and tt production), selected in a region where resummation effects in the coefficient functions are expected to be negligible, since for these observables resummation is not yet available in HELL. Another family is then obtained by including only the DIS datasets in the fit, such that resummation is consistently included for all datapoints. Three variants of these DIS-only fits have been created by enlarging the dataset to include pseudo-data from possible future DIS experiments, namely the Large Hadron-electron Collider (LHeC), the Future Circular electron-hadron Collider (FCC-eh), and both.
For each family, four fits have been performed, with NLO, NLO+NLL, NNLO and NNLO+NLL theory (except for the LHeC and FCC families where only NNLO and NNLO+NLL is available). In all cases, the resummation makes use of the LL anomalous dimension, which, as suggested in Ref. [16] and confirmed in this work, is not the best choice. The new version of HELL released with this work, 3.0, uses the NLL anomalous dimension rather than the LL one as default, so future fits including high-energy resummation will be performed with this new setup. So far, only a single PDF set with resummation at NNLO+NLL has been determined using the NLL anomalous dimension, which has been used in Ref. [1] to investigate the effect of subleading logarithms. However, this set is based on the DIS-only dataset, and as such it suffers from larger uncertainties and it is then not suitable for phenomenology. Nevertheless, its existence will be helpful to investigate the effect of computing consistently the Higgs cross section with our favourite choice of NLL anomalous dimension.
We have to warn the Reader that the HELL 2.0 version of the code [4] used for the aforementioned PDF fits was based on an incorrect resummation formula, which produced spurious NLL contributions to the P gg splitting function beyond O(α 3 s ), and affected other splitting functions and coefficient functions beyond their logarithmic accuracy. The issue has been corrected in HELL 3.0 [16]. The effect of the correction at the level of splitting functions and coefficient functions appears to be reasonably small [16], especially in the kinematic region of HERA, so we expect that the resulting PDFs are not severely affected by the issue. We stress however that the difference between the LL and NLL formulations of the resummation matched to NLO or NNLO is significantly reduced after the correction: therefore, the non-negligible difference in the PDFs [1] obtained with these two formulations using the previous version of the code will likely be reduced significantly in future PDF fits based on HELL 3.0.
In the rest of this section we will proceed as follows. First, we take the global PDF sets of Ref. [1] and compute predictions for the resummed Higgs cross section. Then, we will use the DISonly PDFs to study the impact of subleading terms, both at the level of PDFs and of the coefficient functions. Additionally, in the context of the DIS-only sets we will investigate the reduction of the PDF uncertainty on the Higgs cross section that could be achieved with future DIS experiments.
Let us start with the PDFs based on the global dataset. We consider m H = 125 GeV (physical Higgs) and m t = 173 GeV, and compute the cross section as a function of the collider center-of-mass energy √ s. We set the scales to µ R = µ F = m H /2, which is our default central choice. In Fig. 5 (left plot) we show the cross section at N 3 LO and N 3 LO+LL for a range of collider energies which spans from a Tevatron-like 17 energy of √ s = 2 TeV to a FCC-hh energy of √ s = 100 TeV. Since the cross section changes significantly over this large range of energies, we present the results as ratios (Kfactors) with respect to the fixed-order N 3 LO prediction. For the fixed order (green) and resummed (red) predictions we use the NNLO and NNLO+NLL global PDF sets of Ref. [1], respectively. The uncertainty band shown is due to the PDF uncertainty only. We see that the effect of resummation is small and compatible within the PDF uncertainty for small collider energies, up to the current ggH production cross section ---effect of small-x resummation N 3 LO using NNLO PDFs N 3 LO using NLO PDFs N 3 LO+LL using NNLO+NLL PDFs N 3 LO+LL using NLO+NLL PDFs This huge effect may seem surprising, and thus deserves a careful investigation. First of all, we note that basically the whole effect comes from the use of resummed PDFs, while the effect of the resummation in the coefficient function is almost negligible. Indeed, in the same plot there is an additional curve (dashed blue) obtained by computing the fixed-order N 3 LO cross section with the resummed PDFs: this curve, which differs from the red one only by the resummed contributions to the coefficient function, is basically identical to it, except for a tiny deviation visible only at large collider energies grater than √ s ∼ 30 TeV. We argue that the origin of this huge difference between the predictions obtained with either the NNLO or the NNLO+NLL PDFs is due to the former being unreliable at small x, due to a perturbative instability in the splitting functions at NNLO, in turn due to the unresummed small-x logarithms. Indeed, in Ref. [1] it was observed that the behaviour of the NNLO gluon PDF at small x is rather different from that of the NLO PDF, which in turn is quite similar to both the NLO+NLL and NNLO+NLL resummed gluon PDFs. Namely, the perturbative progression of the PDFs is perturbatively stable at small x when resummation is included, but unstable when resummation is not included, the instability starting to appear at NNLO.
To understand how much of this PDF behaviour is reflected on the Higgs cross section, we show in Fig. 5 (right plot) the fixed-order and resummed cross sections using NLO PDFs (dashed green), NNLO PDFs (solid green), NLO+NLL PDFs (dashed red) and NNLO+NLL PDFs (solid red). We observe indeed that at high collider energies (which probe smaller x and are thus more sensitive to small-x logarithms and their resummation) all curves except the one with NNLO PDFs are grouped together, indicating that the small-x instability of the NNLO is really the culprit of the huge difference between fixed-order and resummed results at high collider energies. Indeed, the resummed result with NNLO+NLL PDFs is a reasonably small correction to the results obtained with either NLO or NLO+NLL PDFs at high energies. We conclude that the effect of small-x resummation on the Higgs cross section is per se not surprisingly large; however, using NNLO PDFs gives rise to unreliable results at high energies, due to the instability at small-x, which is not even covered by the PDF uncertainty. This effect would be even more marked with N 3 LO PDFs, since N 3 LO splitting functions suffer from stronger instabilities, as demonstrated in Ref. [16]. Thus, contrary to the common lore, using N 3 LO PDFs for a N 3 LO cross section such as the Higgs  cross section would produce results which are even less reliable than those with lower order PDFs. Therefore, at high energies precise and reliable predictions can only be based on small-x resummed PDF sets.
We now return to the observation that the resummation in the coefficient function has a tiny effect. This fact is partly due to the fact that we are adding resummation on top of the already rather precise N 3 LO prediction, and is in perfect agreement with the parton-level behaviour observed in Sect. 3.2, together with the observation that the Higgs cross section is threshold dominated. However, we have also noted that the partonic behaviour is rather different when N 3 LO is supplemented with the resummation computed with the LL anomalous dimension. In that case, the resummation has an effect also at larger z. While we believe that this effect is spurious, it is interesting to see how it affects the physical cross section. This is also interesting because the PDF sets of Ref. [1] have been obtained using resummation based on the LL anomalous dimension. 18 (However, the instability of the LL anomalous dimension appears when expanded to O(α 4 s ), which is not the case for the NNLO+NLL resummation used in the PDFs, for which the use of the LL formulation can be considered reliable.) Thus, in the left plot of Fig. 5 we also show the LL version of the resummed prediction (dot-dashed red). In this case, the effect is rather large, even for small collider energies where we expect resummation to have no effect: this is entirely due to the sizeable contribution of the resummation at z ∼ 10 −1 (see Fig. 3), and confirms the spurious nature of such effect. Interestingly, the effect is positive, i.e. it does not compensate in any way the effect of the resummation in the PDFs, which would be expected if the process were included in the PDF fit.
While we have found strong motivations to consider this effect a spurious one, we have suggested in Sect. 3.2 to use the difference of the resummed predictions obtained with NLL and LL anomalous dimension as an uncertainty due to unknown subleading logarithmic contributions. This choice is certainly conservative if one considers the effect on the coefficient function alone, as the uncertainty is much larger than the central effect itself. However, we shall not forget that subleading logarithmic contributions may have sizeable effects in the PDFs as well, which are probably not taken into account by the PDF uncertainty (see also discussion in Ref. [1]). Thus, this uncertainty has the role to also account for subleading logarithms in PDFs, e.g. to compensate for the fact that these PDFs have not been obtained using the NLL anomalous dimension.
It would be interesting to quantify how large the uncertainty from subleading logarithmic contributions in the PDFs can be. One way to do so is to use a PDF set which has been determined using resummation implemented through the NLL anomalous dimension, which clearly differs by subleading contributions with respect to the ones based on the LL anomalous dimension. As anticipated, in Ref. [1] a single NNLO+NLL fit with this setup 19 has been performed, however just in the context of the DIS-only dataset. Thus, to investigate these effects in a consistent manner, we need to consider the DIS-only fits, which however suffer from larger uncertainties and are thus not suitable for phenomenological applications. In Fig. 6 we show (left plot) the resummed cross section (normalized to the N 3 LO one computed with NNLO PDFs) with four different combinations of choices of subleading contributions: using consistently the LL anomalous dimension in both PDFs and coefficient functions (dot-dashed red), using the LL anomalous dimension in the PDFs and the NLL one in the coefficient functions (solid red, our default), using consistently the NLL anomalous dimension in both PDFs and coefficient functions (solid blue), and using the NLL anomalous dimension in the PDFs and the LL one in the coefficient functions (dot-dashed blue). We restrict our attention to the LHC-FCC energy range, and show on our default prediction (solid red) the PDF uncertainty band (darker red area) and the sum in quadrature of it with the "subleading logarithmic uncertainty" as defined above, namely by the difference between solid and dot-dashed red (lighter red area). The solid blue curve is what we would consider the new default prediction, as it uses consistently the NLL anomalous dimension, as suggested in Ref. [16] and here. We note that such prediction is smaller than our default one, reaching "just" a 6% increase over the N 3 LO at FCC, and suggesting that our current default prediction may overestimate the real effect. Nevertheless, we see that our full uncertainty band reasonably takes into account the difference between the two predictions, even though the blue curve lies outside the band for √ s 30 TeV. However, we need to keep in mind that these PDFs are based on the previous version of HELL, where the difference between LL and NLL formulations was larger than in the new corrected version, as we commented before. We may realistically expect that with the new version of the code the PDF sets corresponding to the two variants of the resummaiton be closer to each other, such that our uncertainty band successfully covers such effect. A definitive answer can only be obtained in future, when a (possibly global) PDF fit will be performed with the new HELL 3.0 default, ideally also including the resummation of hadron-hadron collider observables, most importantly Drell-Yan cross sections, which can directly constrain small-x PDFs at larger Q 2 and then reduce an unavoidable source of uncertainty coming from the large portion of DGLAP evolution from the low-Q 2 HERA region to the ElectroWeak scale.
In the context of DIS-only fits, in Ref. [1] it has been studied the impact of the inclusion of pseudo-data from possible future DIS experiments at LHeC and FCC-eh. Both datasets provide a significant reduction on the PDF uncertainty at small x. Most of the reduction is provided by the FCC-eh dataset, with the LHeC dataset providing only an extra little reduction. Being realistic (it is unlikely that both facilities will be built) and also wanting to maximize the effect of the new data, we decide to consider the PDFs obtained with the addition of the FCC-eh dataset alone. In Fig. 6 (right plot) we show the relative PDF uncertainty of the N 3 LO+LL result for the real DIS-only fit (dashed blue) and the futuristic DIS-only fit including FCC-eh pseudo-data (dot-dashed green). We see that indeed the reduction is significant and important in the highenergy region. However, it is way less dramatic than the analogous reduction visible in the gluon PDF (see Ref. [1]). This is due to the fact that we are considering an inclusive cross section, which, according to Eq. (3.19), contains contributions from all the regions of x from τ to 1. Thus, the strong uncertainty reduction on the gluon at small x has only a limited benefit on the full PDF uncertainty of the cross section even at rather large collider energies. Indeed, for comparison, in the plot the 19 However, as already mentioned, also this fit was performed prior to the correction in the resummation code.  uncertainty obtained with the global dataset (and thus without FCC-eh) is also shown (solid red): this uncertainty is always smaller than the DIS-only with FCC-eh one, up to the FCC-hh energy where they become comparable. Thus, for the inclusive cross section, future DIS experiments may lead to an increased precision at high energies, but also precise hadron collider data can. However, when considering differential observables, which are more directly sensitive to the PDFs at specific values of the momentum fraction, the uncertainty reduction provided by FCC-eh or LHeC may be more substantial.
So far we have presented results (with and without resummation) at N 3 LO. To conclude the section, we present some representative results at previous orders. In Fig. 7  At resummed level, the PDF uncertainty (blue), its sum in quadrature with the uncertainty from subleading logarithms (salmon), and their sum in quadrature with the scale uncertainty (yellow). We observe that since most of the effect of the resummation is due to the PDFs, the increase in the cross section is more or less independent of the perturbative order. Therefore the perturbative progression does not improve significantly when adding resummation, 20 even though a marginal improvement is anyway visible -for instance, at the FCC-hh the NLO full uncertainty band does not cover the NNLO result, while the NLO+LL band does cover the central NNLO+LL result. It is interesting to note that the uncertainty from subleading logarithmic contributions is negligible at NLO+LL, small at NNLO+LL and quite large (comparable with scale uncertainty for HE-LHC and FCC-hh) at N 3 LO+LL. Because we compute this uncertainty as the difference between using NLL and LL anomalous dimensions in the resummation of coefficient functions, this pattern shows that these two approaches give quantitatively similar results at NLO+LL and NNLO+LL, but as we have already noted they differ significantly at N 3 LO+LL, in agreement with the parton level results presented in Sect. 3.2. We do not report explicit numerical results, as these have been already presented in Ref. [17], where the contribution from threshold resummation is also included, which is known to stabilize the perturbative expansion of the Higgs cross section, and additional corrections due to e.g. the bottom and charm quark running in the loop are considered.

Conclusions
In this work we have extended the resummation formalism for partonic coefficient functions originally developed for deep inelastic scattering [3] to the case of two hadron in the initial state, relevant for LHC. In particular, at the leading logarithmic accuracy we considered, only processes which are initiated by two gluons at LO, such as Higgs production in gluon fusion, top-pair production, jet production, etc., are non-trivial, while processes which are quark initiated like Drell-Yan resum only through a single initial state leg at this order, and are thus treated identically to the single-hadron case. We have demonstrated the equivalence of our (more general) approach with the original ABF approach of Ref. [9] under specific assumptions, and provided all the ingredients needed to match resummed results to fixed-order computations up to N 3 LO. This formalism has been implemented in the new version of the public code HELL 3.0.
We then studied a specific hadron-hadron collider process, namely Higgs production in gluon fusion. The partonic coefficient functions with incoming off-shell gluons needed for obtaining the resummed on-shell coefficient functions for this process have been computed a while ago [61][62][63]. However, it was possible to obtain consistent resummed predictions only thanks to two recent developments. On the one hand, the creation of the public code HELL which implements the formalism for resummation developed in Refs. [3,4] and extended to the hadron-hadron collider case in this work. On the other hand, the existence of PDF sets which have been obtained using small-x resummation (from HELL) in their determination and evolution [1,2].
Comparing the Higgs cross section at N 3 LO supplemented by small-x resummation using resummed NNLO+NLL PDFs with the (current standard according to the LHC HXSWG) fixed-order N 3 LO prediction using NNLO PDFs, we have found that the cross section increases mildly (+1%) at current LHC energy, and increases more substantially for larger collider energies, reaching +4% at HE-LHC ( √ s = 27 TeV) and +10% at FCC-hh ( √ s = 100 TeV). Almost all of this effect comes from the use of resummed PDFs, and in particular it is due to the fact that NNLO PDFs are unstable at small-x due to the presence in the three-loop splitting functions of large unresummed logarithms of x [1]. The effect would be potentially much larger if (yet unavailable) N 3 LO PDFs were used, since four-loop splitting functions are even more unstable due to larger powers of the logarithms at small-x [16]. The main conclusion that we draw is that predictions based on NNLO PDFs and in future on N 3 LO PDFs will be unreliable for processes which are sensitive to small-x PDFs, due to the bias induced by the perturbative instability of the splitting functions. While for the inclusive Higgs cross section this seems to be the case only at future colliders, for differential observables which are more directly sensitive to PDFs at a given momentum fraction this conclusion will hold also at the LHC in specific kinematic configurations (e.g., large rapidities). In these cases, only predictions based on small-x resummed computation can be considered as reliable.
At the moment, the main limitation of small-x resummation is its limited logarithmic accuracy. For DGLAP evolution, resummation is known at NLL, while for the coefficient functions it is known only at LL. In this work we have also studied the potential effect of subleading logarithmic contributions to the Higgs cross section, by computing different theoretical predictions which differ by subleading terms both in the coefficient functions and in the PDFs. The effect is potentially large, and while the qualitative conclusions of this study remain unchanged, achieving high precision requires the extension of the small-x resummation formalism to higher logarithmic order. This ambitious goal is left to future work.
The new 3.0 version of HELL which contains all these new developments is publicly available for download at the address www.ge.infn.it/∼bonvini/hell It also uses a new default for the implementation of the resummation, as discusses in Ref. [16]. HELL 3.0 has been used in Ref. [17] to obtain double-resummed predictions at threshold (large x) and at high energy (small x) for the Higgs cross section at LHC and beyond.
Thus, the integral over virtualities of a function F (ξ 1 , ξ 2 ) transforms as where in the last line we have assumed F to be symmetric, so that the problematic region ξ 1 = ξ 2 lies at the boundary of the integration domain and can be better integrated numerically. In terms of these variables, the form factors [62,63] have a simpler form given by and . (A.14) All these expressions have been coded in HELL 3.0. In some particular limits, where some of the functions fail to evaluate numerically (mostly due to the square root terms), Taylor expansions are used to overcome this problem.
In the actual definition of the resummed coefficient functions, Eqs. (2.19a) and (2.19b), the integration extends from the position of the Landau pole ξ 0 to infinity, and the off-shell coefficient appears with derivatives with respect to ξ 1 and ξ 2 . The second fact is per se not a problem, except that these derivatives must be computed analytically both for speed reasons and to avoid proliferation of numerical errors. Therefore, it is useful to limit as much as possible the number of derivatives to be computed. To do so, we first observe that we do not need to treat identically the contributions from f 1 and f 2 , Eq. (A.1). Indeed, in our numerical implementation we use the expression in which the derivatives act on the coefficient function for the f 1 contribution, while we use the one with derivatives on the evolution functions for the f 2 contribution. Making the notation very schematic and omitting all arguments except the virtualities, we write Eq. (2.19a) as where U is a shorthand for U ht ABF , and U (ξ) is its the derivative with respect to ξ. The term proportional to f 2 is then treated as described above, namely by changing variables according to Eq. (A.5) and using the symmetry to integrate only for positive y's. The contribution to C gg from f 1 , which we call C 1 for simplicity, is instead manipulated as follows where we have defined U ± ≡ U (t(1 ± y)), (A. 17) and U ± are still derivatives with respect to the full argument. In the first step in Eq. (A.16) we have simply performed the change of variables; in the second step we have integrated by parts some contributions to remove double t and double y derivatives (all boundary terms vanish). At this point one can use the symmetry y → −y to restrict the integration to positive y's up to an overall factor of 2. Eq. (A. 16) is what we use in the code, and provides a stable numerical evaluation of the integral, with the advantage of depending on a single second derivative of the off-shell coefficient function. The auxiliary function Eq. (2.19b) is instead much simpler to treat. First, the f 2 term proportional to |A 3 | 2 does not contribute, since it is multiplied by ξ 1 ξ 2 and one of them is zero (say, ξ 2 = 0), so we have Second, there is a single derivative, which can be directly obtained from ∂f 1 /∂t used above. In fact, the form factor becomes much simpler in the limit ξ 2 = 0, i.e. y = 1, with C 0 (t, 1) = 1 1 + 2t Li being now T + = 1 + 2/(ρ t t). This analytical expression is also useful for cross-checking numerically part of the results used above in the C gg case.
We can now discuss the implication of restricting the integration to ξ 1,2 > ξ 0 = exp −1 αsβ0 . Let us start with the one-dimensional case, C aux , Eq. (A.18). The integrand is peaked at ξ ∼ 1, and drops at large ξ as a negative power of ξ (in this case, as 1/ξ 3 ). Thus, we do not loose precision if we approximate the integrand as where m > 0 cuts off the large ξ region which gives a negligible contribution to the integral, and F (ξ) is a generic name for the integrand. In practice, we have noticed that m = 3 is sufficiently small to guarantee numerical stability and at the same time sufficiently large to keep the important region of the integral and cut away only negligible corrections. We then split the integrand in two pieces, from ξ 0 to 1 and from 1 to ξ −m 0 , and perform the change of variables ξ = exp(−u) and ξ = exp(mu) respectively: such that the integration region is in the unit hypercube (of dimension 1 in this case), and thus directly usable in standard numerical integration routines. This expression is what is used in HELL for the one-dimensional case.
In the two-dimensional case, we need to convert the two conditions ξ 1 > ξ 0 and ξ 2 > ξ 0 into a condition for the y, t integration range. Assuming to integrate in y first, the condition becomes The integral of a generic function F (t, y), once the y integration is symmetrized, can be treated as in the one-dimensional case, approximating the integral and performing subsequent changes of variables, where in the last line we first changed variable y = (1 − ξ 0 /e −u )w in the first y integral and y = (1 − ξ 0 /e mu )w in the second y integral, and then we used again u = −v log ξ 0 . As before, the final result is integrated in the unit hypercube (of dimension 2 in this case), and thus immediately usable for numerical integration. This expression is implemented in HELL for the two-dimensional case.

A.2 Impact factor and its expansion coefficients
We now move to the computation of the coefficients of the M 1,2 expansion of the Mellin transform of the off-shell coefficient function, Eq. (2.20). Such Mellin transform is equivalent to Eq. (A.15) after replacing and letting ξ 0 → 0. We can thus start from Eq. (A. 16) and, after integrating by parts in t the last term, we arrive at (again omitting all non-crucial arguments) 22 Our goal is now to expand this expression in powers of M 1 and M 2 , to construct the coefficients C kj , Eq. (2.21). We observe however that there is a term in Eq. (A.27) which seems to give rise to negative powers of M 1,2 , namely the one with M 1 + M 2 in the denominator. When expanding (tQ 2 /µ 2 F ) M1+M2 in powers of M 1 + M 2 all terms except the zero-th order term will compensate the denominator. Thus, the only term which is potentially dangerous is the zero-th order one, which reads where the coefficients c a,b of the M ± expansion are given by Once these coefficients are known, they can be converted to the desired coefficientsC kj , Eq. suffers from a definition of the coefficientsC kj in terms of integrals that are not easy to perform numerically and give rise to large numerical errors. Therefore, our construction, despite being somewhat involved, has the big advantage of reducing the numerical error significantly, which was possible by exploiting the symmetry of the off-shell coefficient function. We add that the construction presented in this subsection was actually already used for computing these coefficients for Ref. [28], but it is presented in this detail here for the first time.