Differential heavy quark pair production at small x

We consider the production of a heavy quark pair in proton–proton collisions. For bottom and charm quarks, the final state invariant mass is typically much smaller than the collider energy (e.g. at the LHC), so that high-energy logarithms may spoil the perturbativity of the theoretical prediction at fixed order. The resummation of these logarithms to all orders is thus needed to obtain reliable predictions. In this work, we extend previous results on high-energy (or small-x) resummation to differential distributions in rapidity, transverse momentum and invariant mass, and implement them in the public code HELL.


Introduction
In the era of LHC precision physics, considerable efforts are required to match theoretical prediction with experimental accuracy.Such an endeavour requires several different inputs, e.g.high-order predictions for partonic processes, high-quality parton distributions and all-order resummation of large logarithmic contributions.
In this work, we focus on the latter and specifically on the so-called high-energy logarithms of the form α n s 1 x log k 1 x , k < n, where x is a dimensionless scaling variable that becomes small when the collider energy s is large.These perturbative terms arise beyond the leading order in both the partonic cross sections and the DGLAP splitting functions governing PDF evolution (in MSlike schemes).At the energy scales of many LHC processes, x 1 and these logarithms spoil the perturbativity of the fixed-order results.This calls for an all-order resummation of these corrections.
One of the key steps to achieve a consistent resummed prediction is the resummation of partonic cross sections, which can be carried out to leading logarithmic (LL) precision using the k t factorization theorem [26][27][28].Recently, the resummation technique for partonic cross sections has been reformulated and adapted for stable numerical implementation [29][30][31].This led to the release of the High-Energy Large Logarithms (HELL) public code, which aims to provide a systematic framework for implementing small-x resummation.
So far, only inclusive observables have been considered in HELL.The sensitivity of inclusive observables to resummation effects is, however, limited.Indeed, the small-x region at parton level is mixed with the medium-and high-x regions in the convolution that defines the hadronlevel cross section, thereby smoothening out much of the impact of high-energy logarithms (see e.g.Ref. [31]).Differential distributions, instead, can be more directly sensitive to specific values of partonic x, thereby enhancing the effect of small-x resummation in some kinematic regions.Moreover differential distributions are of greater phenomenological interest, as they can be compared more directly with experimental measurements.
In this work we will focus on invariant mass, rapidity and transverse momentum distributions.The resummation of small-x logarithms in these differential cross sections was developed in Refs.[32][33][34], focussing on Higgs production via gluon fusion.Here, we revisit these results and extend them to the modern resummation formalism of Refs.[29][30][31], thereby allowing for a stable numerical implementation thus opening the door to phenomenological studies.
We apply our findings to heavy flavour pair production, and construct resummed predictions for distribution in invariant mass, rapidity and transverse momentum of either the heavy-quark pair or one of the heavy quarks.This process is particularly interesting due to the availability of measurements from the LHCb experiment for the production of charm and bottom quarks in the forward region, where one of the incoming partons is certainly at small x and thus the effect of resummation should be marked.In addition, these data reach values of x down to x ∼ 10 −6 , which is a region of proton momentum fractions so far unexplored, as the HERA data is limited to x 3 • 10 −5 in the perturbative regime.Our results thus provide an important ingredient to refine the determination of PDFs at small-x, which serves both as a test of QCD in extreme regimes and as a tool to improve high-energy phenomenology.All our results are available through the new release of the HELL code.
The structure of this paper is the following.Section 2 is dedicated to presenting the formalism of k t factorization in a proton collider and its use to construct small-x resummed results for differential distributions in the language of HELL.Then, section 3 is dedicated to the application of resummation to differential heavy flavour production, parametrising the final state respectively as the entire quark-antiquark pair or as a single quark.We conclude in section 4, and collect in the appendices various details on analytical expressions for heavy quark production and aspects of numerical implementation.

Multi-differential small-x resummation in HELL
The resummation of small-x logarithms in physical observables is based on k t factorization [26][27][28][35][36][37].The basic observation is that the leading small-x logarithms arise, in a physical gauge, from k t integration over gluon exchanges in the t channel.Therefore, in the small-x limit, the generic amplitude squared can be decomposed into contributions that are two-gluon irreducible (2GI) in the t channel and thus do not contain any logarithmic enhancement.Instead, the small-x logarithms are produced by the integration over the momenta of the gluons connecting these 2GI block.In this way the cross section of the process factorizes [26][27][28] into a process dependent 2GI coefficient, called off-shell coefficient function, and process independent "unintegrated" PDFs that contain the traditional collinear PDFs and the sum over all possible process independent 2GI kernels connected by off-shell gluons.By making explicit the dependence of unintegrated PDF on collinear PDFs and comparing the result with the standard collinear factorization, one finally obtains an expression for the LL resummation of small-x logarithms in the collinear partonic coefficient functions.
The last step of this procedure was traditionally performed in Mellin moment space, which allows to obtain rather simple analytical expressions.Despite the elegance of this result, it was soon realized that subleading effects due to the running of the strong coupling are important and should be included systematically in the resummation procedure to obtain perturbatively stable results [11,38].However, the inclusion of such terms in Mellin space is complicated, and not suitable for efficient numerical implementations.Recently, an alternative but equivalent formulation of the resummation was proposed [29], that solves the technical limitations of the original formulation by working directly in k t space, leading to an efficient numerical implementation.This novel approach is at the core of the public code HELL, and allowed for a number of phenomenological applications [30,31], including the first consistent PDF fits with small-x resummation [23][24][25].
So far, all HELL applications are for inclusive observables (DIS structure functions [29,30] and the total Higgs production cross section [31]).The resummation of differential observables, of obvious interest for LHC phenomenology, has been considered in the Mellin-space formalism.Specifically, resummed expressions for rapidity distributions [32], transverse momentum distributions [33] and double differential distributions in both rapidity and transverse momentum [34] are available.It is the purpose of this section to reformulate these results in the new HELL language, thereby supplementing them with the running coupling contributions and thus providing a ready-to-use numerical implementation.
In this work, we focus on processes at hadron-hadron colliders that are gluon-gluon initiated at lowest order.These include, for instance, Higgs production, jet production, or heavy quark pair production; the latter will be considered as a practical application in section 3. The reason for this choice is that the resummation is simpler, because at LL there are no collinear singularities.In other processes where the lowest order is initiated by (massless) quarks, because small-x logarithms at LL appear from chains of emissions ending with a gluon, the diagram entering the computation of the off-shell coefficient function must contain at least a gluon to (massless) quark splitting, thus producing a collinear singularity.One example is the Drell-Yan process.In such cases, the collinear singularities must be treated at the resummed level (similarly to what is done in DIS, see Ref. [30]).A study of the Drell-Yan process where this issue is addressed at differential level is left to future work [39].
Before moving to the resummation, we establish the notation by presenting the structure of differential distributions in collinear factorization for a process in proton-proton collisions.We consider a generic final-state momentum q (it can be the momentum of a single particle or the sum of momenta of different particles) in the collider center-of-mass frame, and we write the distribution differential in its invariant mass squared Q 2 ≡ q 2 , rapidity Y = 1 2 log q 0 +q 3 q 0 −q 3 and transverse component squared with τ = Q 2 /s (s is the collider energy squared) and the sum extends over all possible partons i, j in each proton.In this expression the function is the parton-level coefficient function, which depends on x = Q 2 /ŝ c (the parton-level analog of τ ) where ŝc is the partonic center-of-mass energy, 1 and on y which is the rapidity of q with respect to the partonic center-of-mass frame, and is related to the proton-level rapidity Y by a longitudinal boost.Indeed, Y − y is the rapidity of the partonic center-of-mass frame with respect to the collider frame, and it is determined by the momentum fractions x1 , x2 of the partons in each proton by Y − y = 1 2 log x1 x2 .Note that we have omitted the dependence of α s and of the coefficient function on the renormalization scale µ R , as such dependence is subleading in the small-x limit we are interested in.Finally, the function is the (collinear) parton luminosity, given by the two PDFs with momentum fractions x1,2 = τ x e ±(Y −y) , and including a θ function which is encodes the condition x1,2 ≤ 1.
Eq. (2.1) can also be rewritten as an integral over the parton momenta x1 , x2 , which represents the direct extension of the analogous formula in DIS.However, this form is more suitable for further manipulations.Indeed, it has the form of a Mellin-Fourier convolution, which implies that it can be diagonalized by taking a Mellin-Fourier transform with respect to τ and Y , where In the last equality we have used the definition Eq. (2.3) and changed variable from x, ȳ to x1,2 = √ xe ±ȳ and used explicitly the θ function to obtain the product of two Mellin transforms We further observe that the dependence on the transverse momentum does not affect the structure of the cross section formula, and thus impacts only the kinematics.

Extension of k t factorization to differential observables in pp collisions
The works of Refs.[32][33][34] provide a proof of a resummation formula for differential observable at LL accuracy with fixed coupling through the so-called ladder-expansion approach.This may seem somewhat different from the original works [26][27][28] where the resummation is obtained by proving a k t factorization and comparing it with the standard collinear factorization formula.In fact, despite the different languages, the two approaches are based exactly on the same underlying factorization property and lead to exactly the same result.It is thus natural to imagine that the results of Refs.[32][33][34] on differential distributions could be reformulated in terms of the k t factorization approach.Indeed, it is not difficult to follow the steps of the derivation of Refs.[32][33][34] and recognise the ingredients of k t factorization to construct a factorized formula.Here, rather than repeating such a derivation, we limit ourselves to formulate the result in k t factorization, showing that it corresponds to the results of Refs.[32][33][34] at LL and fixed coupling.
Similarly to the inclusive case, the differential cross section in k t factorization turns out to be a straightforward extension of the collinear factorization Eq. (2.1) where the partons are replaced by off-shell gluons and integration over this offshellness is added.The result reads dσ dQ2 dY dq and ξ 1,2 = k 2 1,2 /Q 2 are the offshellness of the gluons normalized to the hard scale Q 2 , and k 1,2 are the transverse components of the off-shell gluon momenta (for more details on the kinematics, see App.A).In the expression above dC is the (differential) off-shell coefficient function, representing the process-dependent hard scattering initiated by off-shell gluons.More precisely, it corresponds to the last 2GI part (in the t channel) of the amplitude squared of the process, saturating the off-shell gluon indices with a suitable projector [26][27][28].Everything else is collected into the two unintegrated gluon PDFs F g , that include the standard collinear PDFs and the chain of emissions from the initial parton to the last gluon (the ladder in the language of Refs.[32][33][34]).The integration variables z and η are the analog of x and y of Eq. (2.1), but referred to the center-of-mass frame of the off-shell partons.More precisely, we consider as the parton-level center-of-mass frame in k t -factorization the one obtained if we set the off-shellness equal to zero, so that it is related to the collider frame by a longitudinal boost.More details are given in App. A.
We now show that Eq. (2.7) is equivalent to the result of Ref. [34]. 2 First, we take the Mellin-Fourier transform of this expression with respect to τ and Y , with where we have used the definition Eq. (2.8), changed variable from z, η to x 1,2 = √ ze ±η (the longitudinal proton's momentum fractions carried by each off-shell gluon) and used the θ function to obtain the product of two Mellin transforms At this point we follow Ref. [26][27][28] to write the unintegrated PDF as where γ(N, α s ) is the resummed (gluon) anomalous dimension at LL and R(N, α s ) is a scheme dependent factor.Note that we are ignoring quark contributions for simplicity (we will discuss quarks later in section 2.3).Plugging Eq. (2.12) into Eq.(2.9) we immediately recover the result of Ref. [34].Integrating over q 2 t we also reproduce the result of Ref. [32].To reproduce the result of Ref. [33], which is not differential in rapidity, it is simpler to integrate Eq. (2.7) over Y and then take simply a Mellin transform before using Eq.(2.12).The first step leads to with Because this new rapidity-integrated luminosity has the form of a Mellin convolution, after taking a Mellin transform of the cross section we get Plugging now Eq.(2.12) into Eq.(2.15) we finally obtain the result of Ref. [33].
Because the unintegrated PDF depends on ξ through ξ γ−1 , the integrals over ξ 1,2 take the form of Mellin transforms.Therefore, the results above can be expressed (up to factors) as the γ'th Mellin moments with respect to ξ 1,2 of the partonic off-shell coefficient functions, usually called impact factors.These results can be further supplemented with running coupling effects as described in Refs.[11,38].However, as anticipated, adding running coupling effects to the impact factors is not suitable for numerical implementation.In the next section we will start again from Eq. (2.7) to construct a resummed expression at differential level in the HELL language, which makes the inclusion of running coupling effects straightforward and leads to a stable numerical implementation.

Small-x resummation of differential distributions in the HELL language
The main advantage of the formulation of small-x resummation of Refs.[29][30][31] used in the HELL code is the much simpler and reliable numerical implementation.The reason is twofold.On the one hand, the inclusion of running coupling effects in the resummation can be done straightforwardly without approximation and without affecting the numerical performance, as opposed to the impactfactor approach of Refs.[11,38] where it leads to a divergent series that has to be treated in an approximate way.On the other hand, the result can be expressed in terms of the off-shell coefficient function directly in momentum space, as opposed to the impact-factor formulation where a double Mellin transform in both z and ξ is required for each initial-state off-shell gluon.If these Mellin transforms can be computed analytically, the (very minor) price to pay of the HELL formulation is that the ξ integration has to be performed numerically.However, when the Mellin transform in ξ cannot be computed analytically, the impact-factor formulation becomes problematic, while in the HELL approach this does not represent a problem.
The key step of the HELL approach is to write the unintegrated PDF in terms of the collinear gluon and quark-singlet PDFs in a way that includes running coupling effects.The generic form of such an expression, valid at least at LL, is [29][30][31] where and U (N, k 2 , µ 2 F ) is the evolution function of the collinear gluon 3 from the scale µ 2 F to the scale k 2 , times the scheme dependent function R(N, α s ).In other words, U (N, k 2 , µ 2 F ) is the solution of the DGLAP equation using the small-x LL anomalous dimension, which involves only gluons (they do not mix with the quarks at LL). Keeping running coupling effects when solving the DGLAP evolution equation provides the necessary ingredient to include the sought running coupling effects in the resummation [29][30][31].Conversely, evaluating the evolution function at fixed-coupling we get back Eq.(2.12).
In practice, to simplify the numerical implementation and avoid potential numerical issues, the evolution function is approximated in a way that reproduces exactly the results of Refs.[11,38], namely it is valid at LL and at "leading running coupling" (i.e.leading β 0 terms are retained).Within this approximation it takes the form [30,31] where is a damping function at small ξ, designed to keep unaffected the perturbative expansion of the evolution function while ensuring that it vanishes at the Landau pole ξ 0 as it would do at LL with full running coupling [30], and is the approximated evolution function.The anomalous dimension γ appearing above is in principle the LL anomalous dimension.However, it is convenient to include subleading contributions in it that simply produce subleading effects in the resummation but make the result consistent with the resummation in DGLAP evolution.As this discussion is not central for the present work, we refer the Reader to Ref. [31] for further detail.In the numerical implementation, we will ignore the scheme factor R(N, α s ).The reason is that we use to perform small-x resummation in the so-called Q 0 MS scheme [28,36,40,41] where by definition R(N, α s ) = 1.This scheme differs from the usual MS scheme at relative order α 3 s (at LL), and therefore it can be safely used in conjunction with MS fixed-order computations up to NNLO. 3 It is worth noting that the quark part of (2.17) uses the same evolutor of the gluon part.This is justified as, in the x → 0 limit, the leading splitting functions, Pgg and Pgq, are identical up to a factor C F C A .The subtraction of the δ(ξ) in the quark part, that represents the no-splitting event in which the parton remains collinear, is required as the first splitting of the quark into a gluon must be present, and so that contribution must start at order αs.
Let us focus for simplicity on the gluon contribution only, thus neglecting the quark term in Eq. (2.17).Plugging Eq. (2.17) into Eq.(2.9) we get Comparing this expression with the gluon-gluon channel of the collinear factorization expression Eq. (2.4) and Eq.(2.5) we find the identification So far this is not dissimilar to the approach of older works; in particular, if one replaces U with the LL fixed-coupling expression from Eq. (2.12) one recognises the definition of the impact fator.Here instead, we keep a more generic expression for U and further manipulate the result.Indeed, we notice that the N, b dependence of the right-hand side of Eq. ( 2.23) has the same form of the righthand side of Eq. (2.4) or Eq.(2.9).We thus recognise Eq. (2.23) as the Mellin-Fourier transform of which is expressed as a 4-dimensional integral (to be performed numerically in general) over simple quantities, namely the differential off-shell coefficient function and the evolution factors in physical momentum space.This result is very convenient from a numerical point of view.The two additional integrations over z and η are much simpler to compute than the inverse Mellin-Fourier transform over N and b of Eq. (2.23), especially in HELL, because the anomalous dimension appearing in the definition of U is available in HELL only for values of N along a specific inversion contour, which would not be sufficient here due to the ±ib imaginary shift.Instead, because the evolution function U is universal (process independent), it is computed once and for all in HELL directly in momentum space, and it can be easily used in an expression like Eq. (2.24).Moreover, as already mentioned, with respect to the impact-factor formulation this result easily incorporates the running coupling contributions through the use of the proper evolution function U , Eq. (2.19).
We want to emphasize a difference with respect to previous formulations of resummation in the HELL language.In previous works, because the N dependence of the off-shell coefficient function is subleading, we used to set N = 0 in it before computing the inverse Mellin transform.The main motivation was that the analytical expressions obtained in this way were simpler, and in some cases it is not possible to compute the Mellin transform of the off-shell coefficient function analytically for generic N , but it is possible for N = 0.In our case, this approach would correspond to setting N = 0 in the off-shell coefficient function in Eq. (2.23) before computing the inverse Mellin-Fourier transform.However, when dealing with differential distribution we are often not able to compute analytically the Mellin transform of the off-shell coefficient function, not even in N = 0.So there would be no advantage in setting N = 0 in it.Conversely, there would be disadvantages.Indeed, some physical kinematic constraints would be approximated if computed in N = 0.One of the consequences is that the endpoint of the rapidity distribution, which is a physical property of the process determined by its kinematics, would be wrong when setting N = 0.This is not dissimilar to what has been found in Ref. [30] in the case of DIS, where quark mass effects on kinematic constraints were lost when setting N = 0, requiring a restoration of the constraints by hand.Here, we thus decide that it is much better (and simpler) to keep the subleading N dependence, thereby preserving physical kinematic constraints, without paying any price from the numerical point of view.

All partonic channels
In the resummed expression Eq. (2.24) the key ingredient is (the ξ-derivative of) the evolution function in x space, 4 computed in HELL as the inverse Mellin transform of Eq. (2.19).We observe that such inverse Mellin transform is a distribution.Indeed, expanding U (N, Q 2 ξ, µ 2 F ) in powers of α s the zeroth order term is just 1, whose inverse Mellin is δ(1 − x).Since this is the only distributional contribution in U , we find it more convenient to write it explicitly, ) where U reg is an ordinary function.Computing the ξ-derivative appearing in Eq. (2.17) is not entirely trivial.To do so we first introduce explicitly a factor θ(ξ) in the definition of the evolution function, , which is conceptually harmless as certainly the scale Q 2 ξ = k 2 has to be positive.When deriving we get where in the second step we have used the fact that U reg (N, µ 2 F , µ 2 F ) = 0 and in the last step we have defined the plus distribution according to (2.27) The δ(ξ) term appearing as the derivative of the zeroth order of the evolution has a precise physical meaning: it represents the undisturbed gluon, that does not emit and thus it remains on-shell (ξ = 0).This indeed corresponds to the term subtracted in the quark contribution to the unintegrated PDF, Eq. (2.17).We now observe that the introduction of the plus distribution is not really necessary, because the contribution U reg (N, 0, µ 2 F ) appearing in the first line of Eq. (2.26) is finite.More precisely, because U (N, 0, µ 2 F ) = 0 by construction, Eq. (2.19), we have If this is the case, the first two terms in the first line of Eq. (2.26) would cancel, thus leaving the simpler result which is what we would have obtained if we hadn't introduced the θ function.This implies that the nice physical distinction between the no-emission contribution δ(ξ) and the at-least-one-emission contribution U reg (N, Q 2 ξ, µ 2 F ) gets lost.This is clearly undesirable, and may hint at a problem in the construction of the evolution function.
To understand and overcome this problem, we observe that the ξ → 0 limit of U reg , Eq. (2.28), is localised at large x.But the evolution function at large x is not expected to be accurate, as it is constructed to resum logarithmic contributions at small x.Therefore, we can (and we do) damp the function U reg (x, Q 2 ξ, µ 2 F ) (and thus its ξ-derivative) at large x, with a damping function of the form (1 − x) a (we use a = 2 in the code).After damping, the evolution function satisfies for any value of ξ, including ξ = 0.In this way, we obtain U reg (x, 0, µ 2 F ) = 0 and thus U reg (N, 0, µ 2 F ) = 0, implying that the second term in the first line of Eq. (2.26) vanishes, thus giving In other words, because of the large-x damping, the plus distribution is ineffective.For completeness, we have verified that the numerical integral of Let us now come back to the resummed coefficient function.According to Eq. (2.30), the unintegrated PDF Eq. (2.17) can be rewritten as Physically, the δ(ξ) contribution in the equation above represents the (on-shell) gluon that does not emit, thus producing no logs: this is the fixed-order contribution, and it reproduces the on-shell result.The other term, U reg , is the term containing at least one emission, and thus at least one small-x log.Starting from Eq. (2.31) and proceeding as in the previous section, keeping also the quark contributions this time, we obtain the following expressions These results can be written in a more compact form as and where in the last equation we have used the symmetry ξ 1 ↔ ξ 2 of the off-shell coefficient.So in conclusion the resummed expressions for all channels are written in terms of a "regular" resummed coefficient and two simpler "auxiliary" functions,5 each defined in terms of integrals over ordinary functions (and thus easy to implement numerically).The gg coefficient function also depends on the on-shell limit of the off-shell coefficient; however, whenever the resummed result is matched to a fixed-order computation, this contribution will be subtracted and thus in practical applications it will never be needed.
We observe that the auxiliary functions are obtained by putting on shell one of the incoming gluons.Therefore, they represent a contribution in which resummation, obtained from k t factorization, acts on a single initial state parton, while the other obeys the standard collinear factorization.This resembles the hybrid factorization discussed in Refs.[42][43][44][45][46][47] and used to describe forward production.We believe that our auxiliary contribution does indeed represent the same resummed contributions obtained from the hybrid factorization.However, there may be some differences due to the different approaches to resummation, that we aim at investigating in a future work.

Matching to fixed order
The resummed result Eq. (2.32) contains only the small-x logarithms.For phenomenological applications, it has to be matched with a fixed-order computation.To do this, we need to compute its expansion in powers of α s up to some order, subtract it and replace it with the exact fixed-order result at the same order.
Computing the α s expansion of the resummed result is in principle straightforward, but it needs some care in practice as we shall now see.Note that the α s dependence comes fully from the integrand of Eq. (2.32), and specifically from the evolution function U reg , as the off-shell coefficient function is needed only at the lowest non-trivial order to achieve LL accuracy.To construct the expansion of U reg , let us consider the expansion of U reg first.Because of Eqs.(2.19) and (2.21), it is clear that such an expansion contains powers of log ξ.The first couple of orders take the form (in both the MS and Q 0 MS schemes) having assumed the expansion γ(N, α s ) = α s γ 0 + α 2 s γ 1 + O(α 3 s ) for the resummed anomalous dimension (see Refs. [31,48] for explicit expressions).After computing the derivative with respect to ξ, terms of the form log k ξ/ξ appear.Such terms are not integrable in the ξ → 0 limit, and thus require a regularization procedure.
To do so, we recall that the actual form of the derivative of the evolution function has a plus distribution around U reg , Eq. (2.26).The plus distribution does not play a role at resummed level because to all orders U reg (N, 0, µ 2 F ) = 0, but this is not true order by order.The order by order expansion of the evolution function diverges at ξ = 0, and so the plus distribution becomes essential.
With a slight abuse of notation,6 starting from Eq. (2.36), we can write the first couple of orders of the expansion of U reg , or, in x space, having defined P 00 (x) as the Mellin convolution of two P 0 's, and having used the expansion P (x, α s ) = α s P 0 (x) + α 2 s P 1 (x) + O(α 3 s ) which is the inverse Mellin transform of the resummed anomalous dimension γ(N, α s ).
Plugging these expansions into Eq.(2.34) and Eq.(2.35) we finally obtain the sought perturbative expansion of the resummed result.In particular, we find up to relative x z out of which we can construct the expansion of each coefficient function through Eq. (2.33).
We note in conclusion that this procedure is not dissimilar to what was used in previous works, see e.g.Ref. [31], where the expansion was obtained by expanding the impact factor.However, the derivation obtained here is more "direct", and the result is written in a form that is immediately usable to compute the expansion numerically, without the need to compute analytically the impact factor.

Heavy-quark pair production
Having described the general formalism for the small-x resummation of differential distributions in HELL, we now focus on a specific process: heavy-quark pair production in proton-proton collisions.This process is relevant because at the LHC, and in particular at LHCb, it is measured in the forward region where one parton is at small x, and it can thus provide important constraints on the PDFs (the gluon in particular) in a region of x that is so far unexplored.Moreover, NLO results for this process are available [49,50], and NNLO corrections have also been computed recently [51], making this process suitable for precision studies.
The process can be schematized as where the two incoming protons have light-cone momenta P 1,2 with (P 1 + P 2 ) 2 = s, the outgoing heavy quarks have mass m and momenta p, p with p 2 = p2 = m 2 , and X represents any additional radiation together with the remnants of the protons.For simplicity, we consider the final state to be given by the heavy quarks themselves, thus ignoring their hadronization and eventual decay into lighter hadrons. 7These effects should not affect the impact of resummation, as they factorize (at least at LL) with respect to the hard scattering process.A full phenomenological study of the process including these effects is beyond the scope of this paper and is left to future work.Rather, the scope of this section is to demonstrate the application of the framework introduced in this paper.
The resummation of high-energy logarithms in heavy quark pair production has been considered in the literature, both at the level of the total cross section [27,52] and for some differential observables [45,46,[53][54][55].To perform small-x resummation of differential distributions in our approach, we need to compute the coefficient function of the partonic subprocess where two offshell gluons produce the final state.At lowest order, as appropriate for LL resummation, the process is where the off-shell gluon momenta are parametrized as8 ) In this way, the off-shellness of the gluons is given by a transverse component with respect to the beam axis.The longitudinal momentum fractions x 1,2 correspond to the first argument of the unintegrated PDFs, and their ratio is related to the longitudinal boost of the partonic reference frame used to compute the coefficient function by Note that this frame is not in general the partonic center-of-mass frame due to the presence of a transverse component in the gluon momenta, but it reduces to it in the limit where the gluons are on shell.Additional information on the kinematics is given in Appendix A.
In order to compute the actual off-shell coefficient function, we need to decide which is the vector q with respect to which we want to be differential.There are two natural choices: either q is one of the two heavy quark momenta p or p, or it is the sum of the two momenta, thus representing the momentum of the pair.We now present results for either choice in turn.

Results differential in the single heavy-quark
In this section we consider the final state to be one of the heavy quarks, and thus focus on the differential distribution in the components of the momentum q = p.The details of the computation of the partonic off-shell coefficient function are given in App.A.1.
We start by presenting the resummed result at parton level, computed according to Eqs. (2.33).We consider the resummed coefficient functions for bottom pair production, with m b = 4.6 GeV, double differential in (partonic) rapidity y and transverse momentum p t of the bottom quark. 9In Fig. 1, we show such distribution as a function of y and for fixed p t = 2 GeV, which is a value accessible at LHCb for the production of b-hadrons.In the left panel, we plot separately the regular Eq. (2.34) and auxiliary Eq. (2.35) contributions out of which the various channels can be built according to Eqs. (2.33), while, in the right panel, we combine them according to those equations to construct the coefficient functions for the gg, gq, qg and qq channels.We observe that the shapes of these functions are quite peculiar, mostly due to the peak of the auxiliary contribution at large rapidity.However, we stress that these are parton-level results, and they are not expected to behave smoothly.In fact, due to the all-order nature of these contributions, it is natural that they present some new features missing in the fixed order.
To appreciate the effect of the resummed contributions on physical cross sections, we present the differential distributions after convolution with the PDFs in Fig. 2, considering for definiteness bottom pair production at LHC 13 TeV.We use the NNPDF31sx [23] PDF set that has been obtained in the context of a study on the inclusion of small-x resummation in PDF fits.The advantage of this set is that it provides PDFs consistently obtained with and without the inclusion of small-x resummation.In the following, we will use the same fixed-order PDFs to compute both the fixed-order and the resummed result, in order to emphasise the effect of resummation in the perturbative coefficient.Also, we provide resummed results obtained with the resummed PDFs, to see how much the resummation in PDFs impacts the cross section.However, performing a thorough phenomenological study is beyond the scope of this paper and is left to future work.
The plots of Fig. 2 show the double differential distribution in rapidity Y and transverse momentum p t at various orders (upper plots) and their ratio to the LO (lower plots), as a function of Y and for fixed p t = 2 GeV.In the left plots, we use the same (fixed-order) PDFs for both fixed-order results and resummed results.We show the LO cross section in dashed orange and the NLO one in dashed blue.The latter, obtained from POWHEG-box [56][57][58], is about twice as large as the LO result, which is partly due to the large value of α s at this low scale. 10In solid we plot the LO+LL (orange) and NLO+LL (blue) results.We observe that resummation is a positive correction at LO, of about 50% at central rapidity and decreasing towards the rapidity endpoints.At NLO, the correction of resummation is still positive, but smaller in size, showing that the perturbative expansion converges better when resummation is included.Overall, the NLO+LL result is approximately a 140% correction over the LO across the whole rapidity range except towards the endpoints, where it goes down a bit following the analogous behaviour of the NLO.When the resummed LO+LL and NLO+LL results are computed with resummed PDFs (right plots), the impact of resummation becomes much larger, as a consequence of the fact that the resummed gluon is larger than the fixed-order one at small x [23][24][25].In particular, the NLO+LL curve has a large K-factor at large rapidities, where the contribution from the gluon at small x is dominant.This shows that this observable is very sensitive to the PDFs at small x, and it thus represents an important process to give additional constraints to PDF fits, in agreement with the findings of Refs.[59][60][61]  It is interesting to understand how the various contributions add up to form the resummed result.First of all, we stress that the LO cross section is made of two contributions, one in the gg channel and one in the q q channel.The second one, however, is very small, so the LO curve is almost entirely given by its gg contribution.As far as the resummed result is concerned, we not only distinguish between channels but also between the regular and auxiliary contributions, as given in Eqs.(2.33).The breakdown of the individual resummed contributions to be added to the LO is shown in Fig. 3 (left).We observe that the dominant contributions are those coming from the auxiliary part, both in the gg channel and in the qg + gq channel.The regular contributions are smaller and localised in a region of central rapidity.Also, we note a clear hierarchy in the contributions by the various channels, with the gg dominating over the qg + gq, and the qq being very small.We also stress that the qg + gq channel is symmetric because we plot them together, but each individual contribution, qg and gq, is obviously asymmetric (see Fig. 1).The right plot of Fig. 3 shows the analogous breakdown for the resummed contribution to be added to the NLO to construct the NLO+LL result.The difference here is only in the auxiliary contributions, as the regular contribution starts at relative O(α 2 s ) and is thus unaffected when subtracting the expansion at O(α s ).Because of this subtraction, the auxiliary contributions become comparable with the regular ones at mid rapidities, but they still dominate in the forward region, as expected.
In order to understand the stability of the resummed result, we now discuss its uncertainties.Because our resummed results are accurate at LL only, the first uncertainty we consider is the one coming from the unknown subleading logarithmic contributions.In previous HELL works [ [29][30][31]48] such uncertainty is studied by varying subleading ingredients in the construction of the resummed  anomalous dimension entering the evolution function Eq. (2.19) in two different ways, 11 and by varying the form of the evolution function itself by replacing r(N, α s ), Eq. (2.21) with α s β 0 .The effect of these three independent ways of varying subleading logarithms in the resummed result is then added in quadrature to form a representative uncertainty for the final result.We adopt this procedure here, and we show the resulting uncertainty as a band in Fig. 2. While this way of computing the uncertainties may possibly underestimate the actual size of NLL contributions, it is clear from the plot that the difference between LO+LL and NLO+LL cannot be due to subleading logarithms only, as it is much larger than their uncertainty.Therefore, at this scale and value of p t , contributions that are subleading power at small x are important.This can be seen also by looking at the difference between additive matching (our default) and multiplicative matching, 12shown as a dotted line in the plot.The difference between these two curves, being related to the ratio between the exact NLO and its small-x approximation, also includes the effect of subleading power contributions, and it is indeed outside the uncertainty band from subleading logarithms.
In Fig. 4 we also show the scale uncertainty band of our results.In the left plot, we consider only the factorization scale variation by a canonical factor of 2 up and down, while in the right plot we construct the envelope of the customary 7-point variation of µ F and µ R .Because the rapidity distribution is symmetric, in each plot we show the fixed-order result for negative rapidity and the resummed result for positive rapidity, for a better visualisation of the bands.As far as µ F variation is concerned, we note a clear reduction of the uncertainty after the inclusion of the resummation, demonstrating the perturbative stabilisation that small-x resummation allows to achieve.However, the uncertainty of the resummed result becomes comparable to the one of the fixed order once µ R variations are also taken into account.This is not surprising, for two reasons.The first one is that the value of α s varies significantly as µ R changes because the scale of the process is low (for the same reason, the NLO uncertainty is larger than the LO one).The other reason is that at LL there are no logarithms of µ R in the resummed result to compensate for the change in α s , as µ R dependence in the resummation starts at NLL.To see a reduction of the 7-point uncertainty band 0 1x10  the resummation should be performed at the currently unknown NLL order.
To conclude the section, we now consider the same double differential rapidity distribution but as a function of p t at fixed central rapidity Y = 0.This is shown with fixed-order PDFs in Fig. 5.We observe that going towards large transverse momentum two effects are manifest: the NLO correction grows, and the impact of resummation on the LO gets larger while matching resummation to NLO gives a smaller correction.This suggests that the large NLO contribution at large p t is dominated by small-x logarithms, and once these are resummed the perturbative convergence improves significantly.As the resummation has no direct dependence on the transverse momentum other than in kinematic constraints, this is just a consequence of the kinematics.In particular, we suspect that the smaller available phase space at large p t makes contributions from the low-x region dominant also at central rapidity (at large rapidity this is expected at any p t ).We plan to investigate this effect further in future phenomenological studies.

Results differential in the heavy-quark pair
In this section we consider the final state to be the heavy-quark pair, and so focus on the differential distribution in the components of the momentum q = p + p which is the sum of the momenta of the two heavy quarks.For instance, this choice is appropriate for describing the measurement of a bound state of the heavy quarks, e.g. the J/ψ for cc pairs or the Υ for b b pairs or heavier resonances.The details of the computation of the partonic off-shell coefficient function are given in App.A.2.
Because at the lowest order the process is effectively a 2 to 1 process, the differential coefficient function contains delta functions, Eq. (A.48).This implies that the computation of some of the integrals defining the resummed collinear coefficient functions, as described in section 2.3, can be carried out analytically, partly simplifying the numerical implementation.Explicit expressions are presented in App.B. Note that these simplifications pose some problems in presenting the results.Indeed, for instance, for the triple differential distribution the regular coefficient function Eq. (2.34) is an actual function, while the auxiliary coefficient Eq. (2.35) is a distribution, making a visual comparison at parton level impossible.This problem can be overcome by showing the cross section at hadron level only, after integration with the PDFs.For definiteness, we consider bottom pair production at LHC 13 TeV, with bottom mass m b = 4.6 GeV, as done in the previous section.Similarly, we use the same NNPDF31sx [23] PDF set considered before.
In Fig. 6 we show the triple differential distribution, plotted as a function of the rapidity Y of the pair and at fixed invariant mass Q = 20 GeV and fixed transverse momentum q t = 50 GeV.In this case, the LO curve is not present, as it is proportional to δ(q 2 t ), and so it is zero for any non-zero value of the transverse momentum.Consequently, we cannot show a ratio plot.We observe that the NLO (blue dashed curve) is smaller than the LL curve (solid orange), which is effectively a LO+LL result.After matching with the NLO, the resummed NLO+LL curve (solid blue) represents a small positive correction to the NLO result, pointing toward the still larger LL prediction.This suggests that the inclusion of resummation tends to predict a higher cross section than at NLO, and possibly leads to a better convergence of the perturbative expansion.As we did for Fig. 2, we show on the left the resummed result computed with the same fixed-order PDFs used for the NLO, while we show on the right panel the resummed contribution computed using the resummed PDFs.In this case, the difference between the two options is very mild, probably due to the larger value of τ and the larger invariant mass, showing that this observable is not particularly powerful in constraining the PDFs at small x.Similarly to Fig. 3, we also show the breakdown of the individual contributions to the cross section in Fig. 7.When matching to LO (left plots) we note a pattern similar to what was observed for the single-quark distributions in the previous section.Namely, the auxiliary term dominates over the regular contribution, and the gg channel is larger than the qg + gq, in turn larger than the qq channel.The resummed contribution is positive, consistent with the fact that this pure LL distribution is effectively a LO+LL result.When subtracting the O(α s ) expansion to match the resummation to NLO (right plots) we find a smaller contribution from resummation.Again, we note that the auxiliary contributions are now comparable in size with the regular ones at mid rapidities, but they keep giving a much larger contribution at large rapidities.
We conclude the section by briefly commenting on the uncertainties.In Fig. 6 the resummed curves are supplemented with an uncertainty band computed as discussed in the previous section to estimate the impact of subleading logarithmic contributions.This uncertainty is relatively larger than in the case of single-quark kinematics, but still it cannot account for the full difference between LO+LL and NLO+LL, which thus gets significant contributions from non-small-x effects.The use of multiplicative matching at NLO+LL, probing some subleading power contributions, differs from the additive matching by an amount that is comparable with the uncertainty band from subleading logarithms.Moving to scale variations, we show in Fig. 8 µ F variations on the left plot and a full 7-point variation on the right plot.Considerations similar to what we have done for the singlequark kinematics apply.We limit ourselves to observe that in both plots there is a visible reduction moving from LO+LL to NLO+LL, again hinting at a stabilisation of the perturbative expansion once resummation is included.
Further studies of different distributions and different kinematic configurations, relevant for phenomenological applications, are beyond the scope of this more theoretical paper and are left to future work.

Conclusions
In this paper we have extended the HELL formalism for the small-x resummation of physical observables to differential distributions at LL.We have obtained resummed formulae for differential partonic coefficient functions which are valid for any process that is gluon-gluon initiated at LO.The application of the formalism to other kinds of processes requires the treatment of collinear subtractions to all orders at small x, whose extension at differential level is left to future work [39].
With respect to previous implementations of small-x resummation, we no longer perform an approximation, valid at LL, where the off-shell coefficient function was computed at N = 0 in Mellin space.This approximation simplifies the resummation of inclusive cross sections where such Mellin transform could be computed analytically.Here, because in general we are not able to compute this Mellin transform analytically at differential level, adopting such approximation would not lead to any simplification.Rather, it would break the kinematic limits of the observables, which is clearly undesirable.
We have considered heavy-quark pair production at proton-proton colliders as a representative application of our results.We have resummed distributions differential both in the momentum of a single heavy quark and in the sum of the momenta of both heavy quarks (momentum of the pair).The selection of numerical results presented serves as a demonstration that the methodology works and that it can be used for phenomenology.They also show that the impact of small-x resummation for these observables is significant, as we expect from the low-x values that heavy quark pair production can reach at LHC.However, the results presented in this work do not represent a full phenomenological study, that would also require the description of the hadronisation of the heavy quarks in order to compare with the data.Such a phenomenological study, with the goal to include the process in a PDF fit to improve the PDF quality at low x, is left to future work.
The new version of the HELL code, that implements the resummation of heavy quark pair production at differential level, is available at the url www.roma1.infn.it/∼bonvini/hell The preparation of tables for quick interpolation, needed for phenomenological applications, requires some time and leads to a large amount of data, because of the dependence on many kinematic variables.Therefore, rather than providing general tables within the code (as previously done for DIS and Higgs), we only provide scripts for the generation of such tables, which can then be produced and used directly by the user focussing only on the kinematics of interest.+ 2 where we have used the inverse relations The conditions imposed by the two theta functions in the energies translate easily into conditions on z 1 and z 2 that depend on x 1 and x 2 , namely z From the on-shell conditions Eq. (A.4) we also know that z 1 z 2 x 1 x 2 ≥ 0 and (1−z 1 )(1−z 2 )x 1 x 2 ≥ 0. Because x 1 and x 2 are positive, it follows that z 1 and z 2 satisfy the conditions 0 ≤ z 1,2 ≤ 1, that translate into z(1 After the trivial integration over Q 2 , the phase space can thus be recast as To simplify the notation, we introduce the function and simply write ξ without arguments for short.Putting everything together we have where in the first line 1/2 is the flux factor, σ 0 = 16π 2 α 2 s /Q 2 and the 1/2π comes from the average over ϕ.The matrix element squared |M| 2 is given in Appendix A. 3.
It is most convenient to use the δ function to integrate over ϑ, as all other variables appear at least quadratically.The fact that | cos ϑ| ≤ 1 produces the constraint We then get where the ϕ-differential distribution is evaluated at specific values of ϕ.Note however that integrating over ϕ immediately is not always the best strategy.For instance, when we take the partial on-shell limit ξ 2 → 0 we obtain where the delta functions do not depend on ϕ anymore, thus making its integration trivial This result is needed for the auxiliary function Eq. (2.35), and also for the subtraction of the plus distributions in the perturbative expansion of the resummed result, see sect.2.4.
We can now use these results in the resummation formulae, using the delta functions to perform integrations explicitly when possible.As far as the auxiliary function Eq. (2.35) is concerned, we can start from Eq. (B.5) and use the δ q2 t − ξ 1 to compute the ξ 1 integration.The result is which does not contain any further integration.Note the presence of a delta function in the result, which can be used in the cross section to integrate over parton distributions.The fixed-order expansion of Eq. (B.6) is given according to Eq. (2.40) by We observe that, for this auxiliary function, the expansion is a distribution in q 2 t .This is not an issue: the triple differential distribution is interesting only for non-zero values of q 2 t , and indeed any measurement will require a q 2 t grater than some resolution cutoff.Should one be interested in integrating over q 2 t down to zero, either to compute the integrated distribution or to obtain a binned version of the q 2 t distribution, one has simply to take care of using the plus distribution in the integration.
We now move to the regular function Eq. (2.34).We get where we have used Eq.(B.1) in the first step and Eq.(B.2) in the second step.Note that the theta function inside the integration can be recast in the (physically obvious) constraint The integration over ξ 1,2 shall be performed numerically.The fixed-order expansion of Eq. (B.9) can be computed according to Eq. (2.39).To this end, it is better to start from Eq. (B.8), to obtain having defined This result is not immediately usable as delta functions still appear explicitly, and the cancellation of singularities in ξ 1 , ξ 2 = 0 requires integrating over these delta functions in a proper order.
To do so, we make some observations.First, when ξ 1 or ξ 2 is zero, the coefficient does no longer depend on ϕ in principle.However, when considering the limit ξ 1,2 → 0, a dependence on ϕ remains.So, for later convenience, we keep the last argument having in mind a limit procedure.Second, the last term is proportional to δ(q 2 t ), which is zero everywhere in the q2 t distribution, except for a single point.This point is interesting only for computing the cumulative distribution between q2 t = 0 and some given value, but in this case it is more convenient to consider the integrated distribution and subtract from it the integral from that value to infinity.Therefore, for our purposes, we can assume q2 t > 0 and ignore the last line.Finally, we observe that the integrand is symmetric for the exchange ξ 1 ↔ ξ 2 , as a consequence of the analogous symmetry of the function dC /dϕ.
We can separate the integration region into 4 subregions, divided by the lines ξ 1 = 1 and ξ 2 = 1.As a result we can write 16)  The ξ 1 ↔ ξ 2 symmetry implies I 2 = I 3 , and further allows us to write I 1 as an integral over a triangle.In this way, only one subtraction is needed to make the integral finite; the other one is a finite integrable contribution.As these integrals have to be further integrated in ϕ, one would be tempted to perform this integration first, before proceeding to ξ 1,2 integration.This seems advantageous because the "full" delta function can be easily solved for ϕ (this is what we have already done before) and the subtraction terms are ϕ-independent and thus the integral is trivial.However, proceeding in this way may potentially lead to numerical instabilities.Consider for instance I 2 , where the subtraction is needed to regulate the ξ 1 integral in ξ 1 = 0.After integrating analytically over ϕ, the remaining ξ 2 integration shall be done analytically (using the delta function) for the subtraction term, but numerically (as we have already used the delta function) for the first term.The cancellation between the two terms is then realised after a numerical integration, which may be dangerous.To make the cancellation smoother, it is much safer to use the delta functions to fix the same variable (ξ 2 ) in the first term and in the subtraction term.The same holds for the I 1 integral, in the representation Eq. (B.19).
To do so, we need to solve the delta function for ξ 2 .The zeros of the argument are given by It is then convenient to change integration variable to √ ξ 2 .We get, for a generic function F (ξ 2 ), where the denominator comes from the derivative of the argument of the delta function.According to this result we can find It is easy to check that the explicit integrals in ξ 1 are finite as ξ 1 → 0. Indeed in this limit ξ ± 2 → ± q2 t , so the ξ − 2 contributions die due to the theta functions, and the ξ + 2 contributions become identical to the subtraction terms, thus making the square bracket vanishing in the limit.
The last integral, I 4 , can be performed similarly to I 1 .However, because here there are no subtraction terms, it is possible to use the delta function for any variable, and in this case it may be convenient to do it for ϕ.
Note that I can be simplified as where we have also traded the ξ + 2 > √ ξ 1 condition for a simpler condition on ξ 1 .The second theta function is more stringent than the first one, so the first one can be dropped.The result above can thus be written as  The integral I 4 , having no subtraction in it, can be computed as in Eq. (B.9), using the delta function to fix ϕ.Of course we can use the approach of using the delta function to integrate over ξ 2 also for the resummed result.In this case we find In fact, it is convenient to partition the integration region along the diagonal ξ 1 = ξ 2 , to get

Figure 1 .
Figure 1.The auxiliary Eq. (2.35) and regular Eq. (2.34) functions as a function of partonic rapidity y for single quark production of mass m = 4.6 GeV at pt = 2 GeV and x = 10 −5 (left plot).The resummed coefficient functions at parton level for each partonic channel constructed according to Eq. (2.33) for the same kinematics (right plot).

Figure 2 .
Figure 2. The double differential distribution in rapidity and transverse momentum of the bottom quark, plotted as a function of the rapidity for pt = 2 GeV, for bottom pair production at LHC 13 TeV.The left plots are obtained using NNPDF31sx at fixed order, while in the right plot the resummed result is computed with the resummed PDFs from the same family.The uncertainty band represents an estimate of NLL corrections.

Figure 3 .
Figure3.Breakdown of the individual contributions to the resummed result from the gg, gq + qg and qq channels separating the regular and auxiliary parts.The left plot focuses on the resummed contribution to be matched to the LO, while the right plot focuses on the resummed contribution to be matched to NLO.The results in these plots are obtained using NNPDF31sx with resummation.

Figure 4 .Figure 5 .
Figure 4. Scale uncertainty for the double differential distribution in rapidity and transverse momentum of the bottom quark, plotted as a function of the rapidity for pt = 2 GeV, for bottom pair production at LHC 13 TeV.The left plot shows factorization scale uncertainty only, while the right plot shows the standard 7-point uncertainty envelope.

Figure 6 .
Figure 6.The triple differential distribution in invariant mass, rapidity and transverse momentum of the bottom pair, plotted as a function of the rapidity for Q = 20 GeV and pt = 50 GeV, for bottom pair production at LHC 13 TeV.The left plots are obtained using NNPDF31sx at fixed order, while in the right plot the resummed result is computed with the resummed PDFs from the same family.The uncertainty band represents an estimate for the NLL corrections.

Figure 7 .Figure 8 .
Figure 7. Breakdown of the individual contributions to the resummed triple differential distribution in invariant mass, rapidity and transverse momentum of the bottom pair from the gg, gq + qg and qq channels separating the regular and auxiliary parts.The left plot focuses on the resummed contribution to be matched to the LO, while the right plot focuses on the resummed contribution to be matched to NLO.The results in these plots are obtained using NNPDF31sx with resummation at LHC 13 TeV, as a function of the rapidity, for invariant mass Q = 20 GeV and for transverse momentum qt = 50 GeV.