Threshold Resummation for Hadron Production in the Small-$x$ Region

We study the single hadron inclusive production in the forward rapidity region in proton-nucleus collisions. We find the long-standing negative cross section at next-to-leading-order (NLO) is driven by the large negative threshold logarithmic contributions. We established a factorization theorem for resumming these logarithms with systematically improvable accuracy within the color glass condensate formalism. We demonstrate how the threshold leading logarithmic accuracy can be realized by a suitable scale choice in the NLO results. The NLO spectrums with the threshold logarithms resummed remain positive and impressive agreements with experimental data are observed.

Introduction. Gluon saturation has attracted a lot of attention in recent years in nuclear physics community. This is in particularly true during the rapid development towards the realization of the Electron Ion Collider (EIC), where one of the scientific goals is to search for gluon saturation and to explore the properties of such a regime [1,2]. Gluon saturation plays the key role in understanding high energy proton and heavy nuclei collisions in the high energy limit, where the gluon momentum fraction x is very small. In such a small-x region, the gluon density grows dramatically and enters the nonlinear regime where the gluon recombination becomes equally important to the gluon splitting, and the Color Glass Condensate (CGC) effective field theory [3][4][5] is the proper theoretical framework to describe such a regime. The nonlinear B-JIMWLK equation [6][7][8][9][10][11] replaces the position of the linear BFKL equation [12], which inevitably leads to the gluon saturation [13,14] with a characteristic scale Q s . The saturation scale Q s features the typical transverse momentum of the gluons inside the proton or the nucleus and grows as x decreases.
Experimental efforts have been made to identify the saturation phenomenon.
Earlier experimental hints on gluon saturation include extensive measurements on structure function in deep inelastic scattering at HERA [15], and the strong suppression of both single hadron [16][17][18] and dihadron production [19,21] cross sections at forward rapidity in d+Au collisions at the Relativistic Heavy Ion Collider (RHIC). More recently the measurements at the Large Hadron Collider [22,23] are also compatible with the saturation-model predictions.
In the future, the dedicated measurements at the future EIC will provide further information on gluon saturation.
In order to faithfully and unambiguously establish gluon saturation and its onset, reliable theoretical predictions for the small-x phenomena at colliders are crucial. When Q s Λ QCD and thus the coupling constant α s (Q s ) 1, the theoretical predictions can be built upon perturbative QCD with a suitable factorization frame-work. However, for the semi-hard saturation scale of a few GeVs, α s (Q s ) is typically not small enough. As a consequence, calculations beyond the leading order (LO) are generally required to ensure the convergence of the perturbative results. Recently, tremendous progress have been made in realizing the next-to-leading order (NLO) calculations for the small-x physics [27][28][29][30][31][32][33][34][35].
In the physical processes investigated so far, single inclusive hadron production in proton-nucleus collisions, pA → hX, is among the most studied ones. This will be the main focus of our current paper. The seminal work [29] confirms the CGC factorization for this observable at the NLO. However, the exhibited negative cross section when the hadron transverse momentum p h,⊥ becomes a bit larger was quite a puzzle in the community [36]. Significant efforts have been devoted to understand and resolve this issue, see e.g. [33,[37][38][39][40][41][42] and references therein. In one of the most recent works [33], the approach introduced can maintain the positivity of the cross section to medium p h,⊥ region. However, the cross section eventually becomes negative for even larger p h,⊥ , although such a transverse momentum is perfectly allowed with p h,⊥ √ s. It is thus widely accepted that the practical phenomenological applications of the NLO calculations for this process are by far problematic [43,44].
In this work, we present solid evidence that the threshold logarithm in the QCD perturbation series is the source to the negative cross section. We are able to resum these logarithms to all orders at the leading logarithmic accuracy (LL thr. ). We find that after resummation, the probability of the real soft emission is suppressed and the NLO predictions with the threshold logarithms resummed (NLO + LL thr. ) stay positive and agree well with the experimental data, meanwhile the theoretical uncertainties get dramatically reduced. Early suggestion of such logarithms as solutions to the negative spectrum problem can be found in [42,45]. In the same spirit, it might also be interesting and instructive to notice that collinear logarithms in the NLO BK equation is the main source responsible for the unstable or even negative solutions and an improved equation with these collinear logarithms resummed solves this instability [46][47][48].
Threshold logarithms. Threshold logarithms are common features of the partonic cross sections for hadronic processes [49][50][51]. They are expected to be large and therefore invalidate the truncations in the perturbative expansion in α s , when a massive final state is produced or kinematic constrains are implemented to force the system reaching its maximally allowed energy. Even in cases where all the kinematics are away from the machine threshold, such as the 125 GeV Higgs production at the 13 TeV LHC, the threshold logarithms are still found to be sizable [52], due to the steep falling shape of the parton distribution functions (PDFs) [50], which effectively places a cut-off in the maximally allowed energy and enhances the effects.
The same story happens to pA → hX. The n-th order corrections to the partonic cross section possess the logarithmic structure in the large N c limit where 1 − z = 1 − τ /xξ with x and ξ the momentum fraction in the PDF and the fragmentation function (FF), respectively, as illustrated in fig. 1. Note that 1 − z is the energy fraction carried by the bremsstrahlung radiations.
We have τ = p h,⊥ e y h / √ s, with y h the hadron rapidity and p h,⊥ the transverse momentum. In the forward region, y h is very large.and thus z can quickly approach 1. In this region, the system is reaching the threshold, the radiations can only be soft and the logarithms are large. To make it more specific, we consider the pA → hX at NLO. In the large N c limit, the partonic cross section can be written as [29,33,45,53] where we have factorized out the LO terms. At the same time, c 0 = 2e −γ E with γ E the Euler constant, and p ⊥ = p h,⊥ /ξ is the transverse momentum of the fragmenting parton. We have only written out those (1 − z) singular terms relevant for discussion, but suppress all the (1 − z) non-singular terms for simplicity. Here, X A is the momentum fraction carried by the gluon from the nucleus and X f is the scale due to the rapidity divergence [33,45,55,56].
and W aa is the CGC Wilson line in the adjoint representation. We find it convenient to use the color operator T a i introduced by Catani et al. [54], acting on the i-th parton with color c(c ) in the color space as where T a c,b = if cab if the particle i is a gluon and T a c,b = t a c,b for a final state quark while T a c,b = −t a b,c for a final state antiquark.
When 1−z ∼ O(1), these (1−z) −1 + terms are small and do no harm to the perturbative calculation. In this awayfrom-threshold case, the typical energy scales involved are the longitudinal momentumn · p of the incoming parton moving along n direction where n = (1, 0, 0, 1) and n = (1, 0, 0, −1), and p ⊥ of the out-going parton. The heirachy p ⊥ n · p gives rise to large logarithms lnn ·p p ⊥ , which we will see, can be resumed to all orders by the BK evolution, if the CGC rapidity scale choice X f ∼ X A is made.
However when we increase p h,⊥ , especially in the forward region where y h is large, z quickly approaches its threshold and the threshold logarithmic terms can become extraordinarily large. To demonstrate this point, we plot explicitly this near-threshold situation in fig. 2, using dAu collision at RHIC with √ s = 200 GeV and y h = 2.2 as an example. In the upper panel, the solid curve is the full NLO cross section including the kinematic constraint [29,33,45], while the dashed curve is the full NLO result with the threshold (1 − z) −1 + terms (setting z = 1 in the numerator) in Eq. (2) subtracted. From this comparison, we see clearly that, when the threshold singular terms are absent, the remaining contribution stays positive for the entire p h,⊥ spectrum, while the full NLO prediction quickly drops below zero. In the lower panel of fig. 2, we show the ratio R between the NLO threshold contribution and the full NLO result. To make the plot more evident, we take out the common δ(1 − z) term from both the full NLO and the threshold contributions. We see that for lower p h,⊥ , non-threshold terms are comparable with the threshold contributions. As we increase p h,⊥ , the threshold logarithms soon become overwhelmingly dominant for p h,⊥ > 5 GeV where the NLO cross section starts to become negative, and the ratio R eventually approaches one. Same behaviors are observed in all other forward kinematic settings.
This exercise clearly indicates that 1). the threshold logarithm is the source to the negative cross sections; 2). the threshold logarithm is very large > ∼ 100% of LO in magnitude and thus requires resummation.
Away from threshold. We start with the away-fromthreshold case to introduce our formalism and notations and to highlight how large logarithms are resummed. At LO, the differential cross section within the CGC framework can be written as , with C = N c the number of colors for quark and N 2 c for gluon initial state in large N c limit. We used the LO color space notation |M 0 (b ⊥ ) [54] which includes the CGC (Glauber) Wilson line W icjc (b ⊥ ) with i c and j c the color indices for the in-coming and the out-going partons, fundamental for quark and adjoint for gluon. f i/P is the PDF, x p = p h,⊥ e y h /ξ √ s and D h/j is the FF. Here, ν is the rapidity scale [45,55,56] in our regularization method for the rapidity divergence in the NLO calculations, and will be related later to the gluon rapidity Y A ∼ ln(1/X A ) in the nucleus.
Beyond LO, an all-order factorization theorem can be derived using the machinery of the soft-collinear-effective theory [56][57][58][59] with additional interactions between quarks/gluons and the Wilson line W (x ⊥ ) adding to it [53], which reads Here the collinear function J encodes the corrections from radiations with the momentum scaling as (n · p, n · p, p ⊥ ) ∼ √ s(1, λ 2 , λ), while the soft function S has the momentum scaling k ∼ √ s(λ, λ, λ). The collinear and soft sectors are classified using the observable power counting in [45] and can be calculated perturbatively. At the LO, J (z) = 1δ(1 − z) and S = 1 and we reproduce Eq. (4). Beyond LO, dimensional regularization and additional rapidity regularization are required to regulate the divergences in the collinear and the soft function, which generates the and η poles and the collinear scale µ and the rapidity scale ν dependence [33,45].
With the scale choice µ ∼ p h,⊥ , all logarithms involving the scale µ will be minimized and be absorbed into the usual DGLAP evolutions of PDFs/FFs. Hence we only focus on the logarithms associated with the rapidity scale ν. To all orders, the collinear function J and the soft function S satisfy the rapidity renormalization group equations where F = J or S. The rapidity anomalous dimensions κ γ ν can be read off from the η-poles in the soft and the collinear functions in perturbative calculations. The NLO poles are calculated in [45,53] and lead to with κ = −1(2) for the collinear (soft) function. Here [. . . ] + is the well-known BK evolution kernel, denoted as I BK below. We can solve the renormalization group equation to find F (ν) = U F (ν, ν F ) F (ν F ), and the evolution kernel U F evolves both functions from their natural scale ν F to a common scale ν to evaluate the cross section meanwhile resums large logarithms ln ν ν F . The ν F is determined by minimizing the logarithms in F and leads to ν J =n · p , ν S = p ⊥ for the collinear and the soft sectors [45]. At LL, we find which resums large logarithms of the form ln ν n·p and ln ν p ⊥ in the collinear and the soft functions, respectively.
Here we have used ν/(ν 2 S /ν J ) = ν/(p ⊥ 2 /n · p) = X f /X A , where X f = ν/n · P A and X A = p ⊥ 2 n·pn·P A with P A the momentum of the nucleus, to get the second equation.
The ν-independence of the cross section implies the evolution for the dipole which when traced over, is nothing but the BK equation. Here the color notation T a b,c is defined previously in Eq. (3) and not to be confused with the fundamental representation t a b,c . With the evolution in Eq. (8), the choice of the rapidity scale X f (or equivalently ν) could in principle be arbitrary, since all large logarithms are resummed. One natural choice is to set X f = X A which is nothing but the conventional CGC scale choice. In such a way, one only needs to evolve the CGC dipoles W † ⊗ W since the evolution U J U S = 1 gets eliminated. In other words, all large logarithms ln p ⊥ n·p are effectively absorbed into the dipole evolution, if X f ∼ X A , when away from threshold.
Near threshold. When near the threshold, real energetic collinear radiations are forbidden, since the longitudinal momentum of the emitted gluonn · p(1 − z) → 0 is restricted to be soft as z → 1, while virtual collinear corrections are still allowed [45]. Therefore, in the threshold limit, the collinear function J thr. will only contains the collinear virtual corrections. All real radiations are now soft and encoded in S thr. . In the large N c limit, it is found that still only the soft and collinear modes can contribute at the leading power [45] and we find the form of the factorization theorem remains the same as Eq. (5) but with the replacement J (z) → J thr. and S → S thr. (z).
The NLO threshold collinear function J thr. is exactly the NLO virtual corrections in J , which gives rise to the evolution of the threshold collinear function at LL where I BK,v (r ⊥ ) = e ip ⊥ ·r ⊥ r ⊥ 2 + e −ip ⊥ ·r ⊥ , is the NLO virtual correction to the BK kernel and ν J ∼n · p to avoid the occurrence of the large logarithms within J thr. . On the other hand, the calculation of the NLO threshold soft function is depicted in [45] and is found to be where I BK,r = I BK − I BK,v is the real contribution to the BK evolution kernel. Here S (1) is the NLO soft function for the away-from-threshold case [45,53], which contains the kinematic constraints. The second term got its contribution from the initial and final parton splitting, which will be absorbed into the threshold evolution of the PDFs and the FFs. We can perform the Mellin transformation dzz N −1 S thr. (z) to the soft function to find ν S ∼ p ⊥ ∼n ·p N e −γ E which minimizes the logarithms in S thr. . We find the associated evolution gives We merge both the evolutions to find where we notice that the second term is identical to the away-from-threshold evolution while the additional first term arises to resum the threshold logarithms. The probability for emitting a soft parton (real correction) is suppressed after resummation. The results could also be realized by considering strongly rapidity ordered soft emissions [53] and are extensible to other small-x processes. From the result, we see that, when near threshold, suppose we still stick to the conventional CGC scale choice X f = X A , then there requires an additional evolution factor to take care of the threshold impacts which can not be covered by simply evolving the CGC dipole.
Instead, we can dynamically determine X f by demanding it minimizing the exponent in Eq. (12) following the similar procedure in [60,61], and hence eliminate the complicated evolution but still maintain the threshold resummation to all orders. The idea is similar to set X f ∼ X A in the away-from-threshold case. We will use this approach for phenomenology studies.
Phenomenology. Now we illustrate the numerical NLO+LL thr. predictions for the kinematics relevant to both the RHIC and LHC experiments. We include all partonic channels. We used MSTW2008 PDF sets [62] for the initial proton and DSS parametrizations [63,64] for the FFs. The CGC dipoles are obtained by solving the LL BK equation with the running coupling correction [65][66][67], with the parameters used in [68]. We set the collinear factorization scale µ = p h,⊥ . For fixed kinematics, we determine the central rapidity scale by scanning through X f (or equivalently ν) numerically to find the value that minimizes the exponent in Eq. (12), which accounts for the threshold LL resummation.
We present the predictions in fig. 3, where we compare the theoretical results with the experimental data in the

FIG. 3. Data versus theory predictions.
forward rapidity region from the charged hadron production in p+Pb collisions at LHC with √ s = 5.02 TeV, 1.5 < y h < 1.8 [23] and the negative hadron productions in d+Au collisions at RHIC with √ s = 200 GeV, y h = 2.2 and y h = 3.2 [16]. From fig. 3, we see that the NLO+LL thr. results (red solid lines) stay positive all the way to large p h,⊥ and show no signs of turning negative. The uncertainty bands represent the rapidity factorization scale X f variation, which are obtained by varying X f around its central value up and down by a factor of 2. We see that the uncertainties are substantially reduced when we go from LO (orange bands) to NLO+LL thr. (red bands). The NLO+LL thr. calculation does an impressive job in describing all the experimental data (green dots in fig. 3). The central values of the NLO+LL thr. predictions slightly overshoot the LHC data for small p h,⊥ but still within errors. The situation is expected to be further improved if a global fit beyond LO is performed to determine the CGC dipole initial condition.
Conclusions. In this paper, we identify the threshold logarithms to be the main source responsible for negative cross section problem of the forward hadron production in proton-nucleus collisions, pA → hX, within the smallx formalism. We develop an all-order factorization theorem with systematically improvable resummation accuracy. We present detailed derivation and numerical study for threshold resummation at LL. We find that the LL thr. resummation can be realized simply by a suitable rapidity scale choice in the NLO calculation. After resummation, all predicted p h,⊥ spectrums are found to be positive all the way to the kinematic boundaries. We compared the NLO+LL thr. results with the LHC and the RHIC data and observed excellent agreements with greatly reduced scale uncertainties, in comparison with the LO results. Our results are ready for more phenomenological applications at the LHC and RHIC, such as global fitting studies of the CGC models beyond LO. Given the universality of the LL thr. structure in hadronic processes, we expect our approach is applicable to many other practical applications of high order CGC predictions for the small-x collider phenomenology.