Threshold Resummation at NNLL for Single-particle Production in Hadronic Collisions

We develop threshold resummation for single-particle inclusive cross sections in hadron-hadron collisions to the level of next-to-next-to-leading logarithm, up to full matching with two-loop hard functions. We define and calculate all one-loop soft functions for all partonic channels. This enables us to separate the hard and soft functions at one loop. Along with these results, the one-loop finite parts of jet functions are used to check that the full soft, collinear and virtual corrections are reproduced to one loop for all partonic reactions. We exhibit these NLO results explicitly. NLO expansions of the resummed cross section match the exact NLO results extremely well numerically, and two loop expansions result in substantial corrections over many kinematic configurations. Explicit results are given in Mellin moment space, and a number of options for generating resummed cross sections are discussed.


Introduction
Single-particle inclusive (1PI) cross sections are among the fundamental processes in QCD, with factorization and evolution properties that are basic results of quantum field theory [1,2]. It is the purpose of this paper to study the resummation of 1PI cross sections for hadrons in hadron-hadron scattering at the level of next-to-next-to-leading logarithmic (NNLL) resummation at partonic threshold [3,4]. We will derive explicit NNLL-resummed partonic hard-scattering functions for all partonic 2 → 2 reactions in terms of moments with respect to the variableŝ 4 , which characterizes the kinematic distance to partonic threshold for the production of a final-state particle with observed transverse momentum and rapidity.
In hadronic scattering at NLO [5,6], large corrections to single-inclusive cross sections are primarily associated with the kinematics of partonic threshold [3,4,7,8,9]. The resummation of such threshold corrections for 1PI cross sections was systematized in Refs. [7,8,10,11,12] for the case of prompt photons, where the resummation was carried out at next-to-leading logarithmic (NLL) level. In particular, the NLO expansion of the NLL resummed cross section can readily be compared to the full NLO result [13,14] for this process. Resummation for single-inclusive hadron production was investigated at NLL for hadronic scattering in Ref. [15] and for photoproduction in Ref. [16]. Especially large corrections from resummation were found in [15] for the case of the rapidity-integrated cross section. Important extensions of resummation techniques for 1PI cross sections in hadronic scattering to next-to-next-to-leading logarithm (NNLL) were made for hadron production in [17] in conventional "direct" QCD, and for prompt photon cross sections in [18] and top production in [19] using soft-collinear effective theory (SCET). The effective theory treatments provide similar results in a Mellin transform space for the partonic cross section, as described below, but differ in their implementation of the transform back to the factorized cross section. Here we will follow a direct QCD approach, although we will make contact with the effective theory results of Ref. [18]. We will rederive the necessary results found in Ref. [17] and the analysis of wide-angle soft radiation in [19], develop further the all-orders factorization properties of the relevant cross sections in terms of direct QCD matrix elements, and provide a formal moment inversion prescription. In addition, we will provide explicit expressions for the relevant hard-scattering and soft functions necessary to describe 1PI hadron production cross sections at NNLL level in hadronic scattering, and exhibit finite NLO contributions that enter the resummed NNLL cross section in moment space with the running coupling evolved to a soft scale.
The path taken in this paper follows the general lines of threshold resummation for dihadron pairs in Ref. [20] at NLL and [21] at NNLL. In the development of the formalism, however, we find several significant differences from the dihadron case. In particular, we conclude that a direct application of the inverse transform to momentum space following the method of [22] is not as appropriate for 1PI cross sections as for the dihadron case. This is manifested by potentially large values of the 1PI cross section in an unphysical region, even though these contributions remain formally power-suppressed in moment space [22,23]. We leave for future work the development of an extension of the method of Ref. [22] to this case. We believe that the studies in this paper will be relevant to classic fixed target data on 1PI cross sections for hadron-hadron scattering, as well as to higher-energy collider data. Extensions of the formalism to photoproduction at NNLL, both resolved and direct, are straightforward.
For very high-energy colliders, jet inclusive cross sections are of special interest, and are clearly related to the results presented here. Fixed-order jet cross sections have been brought to the level of two loops, [24,25], and the impact of threshold resummation has been discussed at least to NLL accuracy in direct QCD [26,27,28,29,30,31,32] and effective theory formalisms [33,34,35,36,37,38]. Inclusive jet cross sections share the underlying kinematics of single-particle inclusive cross sections, and we anticipate that the formalism for single particles described in this paper will have useful applications to single jets.
We begin with a review of the kinematics and the factorization properties of single-particle inclusive cross sections in Sec. 2. The "refactorization" relations upon which threshold resummation for 1PI cross sections is based are reviewed in Sec. 3, where we emphasize similarities and differences compared to dihadron and related cross sections for which threshold resummations have been carried out. In Sec. 4, we rederive results for the relevant jet functions in moment space. The definition of the direct QCD soft function turns out to require a slightly modified treatment of soft gluon phase space, leading to a soft function that differs from that for dihadron cross sections [21], for example. This construction is the subject of Sec. 5. The resulting calculations for the one-loop soft function are given for qq scattering in Sec. 5.3 and summarized for all other partonic processes in an Appendix. The determination of the finite matrices that describe the soft function for all 2 → 2 partonic processes is a main result of this paper. We collect the full resummed cross section in moment space in Sec. 6.1, and in Sec. 6.2 we confirm that all singular behavior at threshold is reproduced at NLO by the one-loop expansions of our hard, soft and jet functions. The results of Sec. 6.1 are presented in an alternate form in Sec. 6.3, with fixed jet and soft evolution scales, confirming that they are consistent with those of [18] when specialized to prompt photon production. We review in Sec. 6.4 several approaches to the application of our results to quantitative cross sections. Section 6.5 presents exploratory numerical tests of the fixed-order (NNLO) expansion of the resummed cross section obtained in our formalism. In an additional Appendix, we provide for completeness the explicit NLO singular contributions at threshold, which are reproduced for all partonic processes by the resummation studied here.

Factorization and kinematics
The classic, collinear-factorized form of the single-particle inclusive cross section pp → hX is with p T the transverse momentum of the observed hadron h, and η its rapidity. The partonic inclusive hard-scattering function ω ab→c , which is accessible to perturbative QCD, describes the production of parton c with transverse momentump T = p T /z c and rapidityη in the partonic center-of-mass frame, the latter related to η bŷ Parton c subsequently fragments into hadron h. We define, for hadronic and partonic kinematics, with s andŝ = x a x b s the hadronic and partonic center-of-mass energies squared, respectively. In Eq. (1), µ F is the factorization scale that ties together the partonic inclusive hard-scattering functions and parton distributions f a,b and fragmentation functions D h c . We omit dependence on the scale at which the perturbative coupling is evaluated, normally denoted by µ R . In principle, ω is independent of this choice, although of course to any fixed order, dependence on µ R appears at the next order.
The conventional partonic variable for 1PI resummation is defined in terms of the momenta p a , p b , p c in the partonic reaction a + b → c + anything aŝ which, as indicated, is the square of the invariant mass of all radiation additional to p c in the partonic final state. In the partonic threshold limit, this quantity vanishes, and the hard-scattering function becomes singular.

Moment analysis of the inclusive hard-scattering function
We write the cross section in Eq. (1) in short as where we have for simplicity omitted all dependence on µ F . Using Eq. (4), the second argument in the inclusive hard-scattering function ω ab→c may also be written aŝ In (5), the lower limit of the integration over the fragmentation variable z c is given by so that z c > z ensures thatŝ 4 > 0. In these terms, we may also write the cross section (5) as so that the last integral takes the form of a genuine convolution of the fragmentation function with the hard-scattering function, at fixedη, which we denote as Moments of Ω ab→c then factor into simple products: As we review below, these are precisely the moments that the resummation formalism gives us near partonic threshold, whereŝ 4 → 0 or alternately N is large in Eq. (11). After resummation, we will perform the Mellin inverse of (10), with C a contour to the right of all singularities ofω. In principle, the moments of the fragmentation functions fall off sufficiently fast so that the integration can be carried out numerically without problem in Eq. (12). So long as the short-distance function has support only for z < 1, as is the case for any finite-order expansion of the resummed expression, the procedure is unique. Assuming that we may use the standard Mellin contour for (12), which, when tilted to the left, acquires an ever larger negative real part along the contour, the Mellin integral will converge very rapidly. The result of this procedure can then be convoluted with the parton distributions in Eq. (5), to give where we have exhibited the ranges of partonic fractional momenta for the incoming hadrons in terms of the quantities x T and η that define the cross section.
We note that the procedure described here for moment inversion is different from the one used for dihadron production considered in Refs. [20] and [21]. For those processes, it was convenient to take a double transform: a Mellin transform in the ratio of the pair invariant mass to the center of mass energy and a Fourier transform in the average rapidity of the observed particles. After resummation, the double inverse transform may be carried out numerically as in Ref. [11]. The reasons for following a different path here will become clear from the refactorization discussion in the following section, and in the explicit form of the moment-space resummation in Sec. 6.

Refactorization in momentum space
Partonic threshold is reached whenx T coshη = 1 and the 2 → 2 partonic subprocess is elastic, leaving no energy available for radiation. This restriction of phase space leads, as usual, to plus distributions in the partonic variable, 1 − y ≡ŝ 4 /ŝ, the most singular of which are double logarithmic, α n s ln 2n−1 (1 − y)/(1 − y) + . Threshold resummation organizes these leading, and in principle all nonleading, distributions of 1 − y at fixed rapidity for the observed hadron.
For such a single-particle inclusive cross section we have a further factorization of the inclusive hard-scattering function ω ab→c near partonic threshold, whereŝ 4 vanishes. This is the lightlike limit for all radiation that accompanies the fragmenting parton c in the final state of the partonic subprocess. As shown in Refs. [7,8], near partonic threshold each particle in the final state may be associated with one of five factoring functions. The first two, labelled here J in , are associated with the two incoming partons. The third, S ab→cr , is a matrix in the space of color tensors, and generates coherent, wide-angle radiation. The fourth, J (r) rec , describes a jet of particles recoiling against the direction of p c in the partonic center of mass frame, and the fifth, J (c) fr , describes radiation collinear to the fragmenting jet itself. The key to threshold resummation in this context, and elsewhere, is that in computations of the inclusive hard-scattering function, ω ab→c , every quantum of final-state radiation can be associated with one of these functions. This feature is directly related to the definition ofŝ 4 in Eq. (4) by where k a and k b are the momenta of particles associated with the incoming jets, k S with the function that generates wide-angle radiation, k r with the jet recoiling against the observed particle, and finally k c with radiation collinear to parton p c . The natural metric for the approach to partonic threshold isŝ 4 , defined in Eq. (4) as the invariant mass associated with Eq. (14).
The next key observation is that in the limit of partonic threshold, the only momentum with fixed, nonvanishing components is k r , associated with the recoil jet. Denoting λ ≡ŝ 4 /ŝ, to leading power in λ,ŝ 4 = (p a + p b − p c ) 2 may be written as a sum, This relation produces a kinematic convolution near partonic threshold between the jet and soft functions into which the cross section factorizes. Of course, there are ambiguities in the definitions of the functions that generate this radiation, especially for wide-angle soft radiation, in the contribution of each function to a portion of phase space. It is precisely this issue that is resolved by the arguments for factorization in Refs. [7,8], for example, and is built into the effective theory treatment of prompt photons in Ref. [18].
We will use the formalism developed in Ref. [7], reexpressed in covariant gauges for threshold resummation in Refs. [20,21]. In this formalism, the incoming jet functions J in develop logarithms not in k r · k a,b directly, but in variables w a and w b , determined by the corresponding terms in (15) through the kinematic relations The fragmenting jet is a function of the fraction of momentum available for collinear radiation, and is denoted by J The recoil jet is a function of its invariant mass directly. We denote it as J (r) Finally, the soft function depends on the ratio with β r a light-like velocity vector in the direction of the recoil jet. Although the variable w S vanishes when k S is collinear to β r , the soft function is constructed so that there are no enhancements in this limit. We can take the all-log resummations for each of these functions from the literature, and we will describe them below.
In the notation we have just reviewed, we can write our re-factorized short distance function as [7,39] where H is a matrix in the space of color tensors, appearing in a trace of possible color factors with the soft matrix S [39]. The matrix H serves as a purely short-distance, non-radiative coefficient function for the soft matrix and jet functions. We will generally refer to it as the "hard function" below, although it should not be confused with the full, inclusive hard-scattering function ω ab→c , of whose refactorization H is a part. Both ω ab→c and H can be computed order by order in perturbation theory. The inclusive hard-scattering function ω ab→c is known explicitly to NLO [5,6], while the hard function H has been computed to NLO in [40,41,21], and even to NNLO in [42].
As indicated, corrections to the refactorized form given in Eq. (20) are not power-divergent at partonic threshold,ŝ 4 = 0, although they may behave logarithmically inŝ 4 /ŝ. The freedom associated with the factorization of ω ab→c in Eq. (20) in terms of four jet functions, a hard matrix, and a soft matrix is encoded into a "refactorization scale", labelled here by µ rf . The hard scattering function, however, is independent of µ rf , For matching comparisons to fixed-order calculations, and other purposes associated with numerical evaluation of the cross section, we will introduce a standard "renormalization scale", which in principle can be chosen independently. For example, to match to an NLO calculation at a renormalization scale like p T , we expand each coupling α s (µ rf ) in terms of α s (p T ). Notice that every function on the right-hand side of Eq. (20) is formally independent of the choice of µ R in these terms, because α s (µ rf ) must take the same value when reexpanded in terms of the coupling at another scale.
An important feature of Eq. (20) is the dependence on the center-of-mass rapidity,η, in the hard and soft functions. Through Eq. (2), this means that these functions depend directly on the observed rapidity of the particle. For dihadron production [21], in contrast, hard and soft functions depend only on boost-invariant differences in the rapidities of the hadrons. As a result, for those processes it is straightforward to perform a Fourier transform in the average rapidity, with η i the physical rapidity of particle i, as mentioned at the end of Sec. 2. In principle, such a transform can be carried out here as well, following the example of Ref. [11] for prompt photons, but the analytic forms would be more challenging for the exponentiated soft functions, which are matrices in this case. Similarly,η appears implicitly in the delta function that defines the convolution in Eq. (20), so that η-dependence enters indirectly in the Mellin moments. These differences led us to forego the use of a Fourier transform in rapidity as in [11], instead using the calculation of the function Ω in Eq. (13).
With the exception of the hard function, each of the factors in Eq. (20) can be computed from universal matrix elements in full QCD, which will be identified below, and they all obey evolution equations that control their behavior in the kinematic limits that provide logarithms ofŝ 4 /ŝ. The hard function H is computed directly from the virtual corrections to the 2 → 2 scattering processes ab → cr, decomposed in terms of color tensors [21,40,41,42]. The refactorization described here has precise analogs in the formalism of soft-collinear effective theory applied at this level for prompt photon production in Ref. [18].
In our analysis we use the feature that convolutions like the right-hand side of Eq. (20) are factorized into products in Mellin-or equivalently Laplace moment space [7,8,17] through for the functions in Eq. (20), each of which has distributions of the form [ln n w/w] + . Corrections to the second equality are suppressed by inverse powers of the moment variable N . Note that the kinematic argument of the soft function is w S √ŝ /µ rf . This dependence is also found, for example, in inclusive electroweak, dijet [26] and dihadron [20,21] cross sections. In this notation, theŝ 4 moments, Eq. (11), of the factorized hard scattering function (20) are now products in terms of hard-scattering and soft functions, and initial state and final state jet functions, identified above, The power of N in each transformed function is determined by whether the corresponding weight w i appears withŝ or √ŝ in the momentum-space refactorized cross section, Eq. (20). On the left-hand side, we have suppressed the argument describing the factorization scale dependence for simplicity. For the initial-state jets the moment variables N i are shifted by simple kinematic factors from the value N in the definition of the moment: which, as noted above, then depend on the center of mass rapidity. Here and below, it is convenient to use the common notation By analogy to the methods employed in Refs. [20] and [21] for dihadron production, we treat the hard-scattering function in Mellin moment space as an intermediate step. Once resummed in Mellin space, we can invert its combination with the fragmentation function to momentum space, as in Eq. (12). To derive the single-particle inclusive cross section, we may then do the resulting integration over x a and x b with parton distributions, as in Eq. (1). We will return to this and other possible prescriptions in Sec. 6.
We next review the resummed forms that we will use for the functions appearing in the moment-space hard scattering function (24), starting with the jet functions. We note that the hard functions H ab→cr in Eq. (24) were already given to one-loop for all partonic subprocesses in Ref. [41] and in our previous paper [21] on dihadron production, and we do not present them again here. The same functions appear in single-inclusive production. Very recently, also the two-loop corrections to the H ab→cr have become available [42], and we will make use of this information in our exploratory phenomenological results.

Evolution in the refactorized cross section
As usual, the functions appearing in refactorized hard scattering cross sections satisfy evolution equations, which control dependence on logarithms of the moment variable. The logarithmic corrections in the moment variable N in each of the functions in the refactorized inclusive hard scattering function, Eq. (24) can be resummed at leading power in N [3,4]. Such resummations have been carried out for various cross sections and in different notations and normalizations, and we shall not try to present a full history here. For our purposes, it is convenient to observe that when the refactorization and factorization scales are taken as equal, each of our jet functions satisfies a specific evolution equation. We will give here the relevant evolution equations for the "in" and "fragmentation" jets (the same equation for all three), for the recoil jet and then the soft function.
First, consider the "in", and "fragmentation" jets, where for consistency of notation we keep an argument 1 = µ F /µ rf . The equation is of a standard form, with γ J (i) the "Drell-Yan (DY)" anomalous dimension [43], When µ F varies at fixed µ rf these same functions satisfy an equation that generates the "cusp" part of collinear evolution, We show how to satisfy this pair of equations in the following section on jet functions. We already see here, however, that the ln N -dependence in the first of these equations is due entirely to the change in the factorization scale.
Next, for the recoil jet, we will show in the following section that this function depends on two ratios involving the moment variable N . The result that we use here is that its dependence on µ = µ rf is independent of N , although the function itself still depends on N , with r the flavor of the parton that initiates the recoil jet. We will find the single-log anomalous dimension γ J (r) rec below. Note that the recoil jet is independent of the factorization scale µ F since it does not involve any parton distributions or fragmentation functions.
The soft function satisfies a familiar equation, with a matrix Γ of anomalous dimensions and without an explicit logarithmic term, [39] The soft function knows nothing about the factorization scale, and the scale µ here is again the refactorization/renormalization scale.
Given the above, we can derive equations for the hard functions H. First, we demand that the short-distance inclusive hard-scattering function be independent of the refactorization scale µ rf at fixed factorization scale, µ F , Second, at fixed µ rf , the function must obey a standard evolution in µ F , which at leading power in N we can write as Here the terms P i,δ are the coefficients of δ(1 − x) in the diagonal evolution kernels for parton i.
The evolution equations for the hard function now follow from the µ rf -independence of ω ab→c , Eq. (31), combined with the evolution of the jet functions and soft matrix, Eqs. (27), (29) and (30). In particular, ln N dependence cancels in the derivative of the hard function H with respect to the refactorization scale at fixed factorization scale, as it must, since the hard scattering function is purely virtual at partonic threshold. This leaves for the µ rf evolution equation for H, where we have dropped the subscript ab → cr labeling the process. To find µ rf -dependence in fixed order computations of H, an equivalent evolution equation is Here we have exhibited the µ rf dependence of the coupling through Meanwhile the dependence of the hard function on the factorization scale is determined by so that That is, we must associate the part of evolution due to parton self-energies with the hard function. Taken together with the purely eikonal evolution equations for the incoming and fragmenting jets, (28), this ensures that the full µ F dependence of the complete perturbative function, ω ab→c , is given by Eq. (32).
The relationship of this kind between factorization and evolution equations was explored in Ref. [44], and equations of precisely this form are a familiar feature of factorized jet and soft functions in direct analyses in QCD [43,45] and in soft-collinear effective theory (SCET) [18,46]. Technical comparisons of the formalisms were given in [45], [47], [48], and [49]. Solutions to Eqs. The solutions to Eqs. (27) -(30) that we will use below relate the jet and soft functions evaluated at a "hard" scale, µ h , of orderŝ to its value at a "soft" scale, µ s , which can, but need not, be chosen to controlN dependence. The hard scale is normally chosen at a value that characterizes the hard scattering. Whatever the choices of µ h and µ s introduced in moment space, it is possible to invert the transform, so long as µ s is large enough that the integral over µ avoids any singularities associated with the running coupling, the "Landau poles". In SCET the conventional choice for the lower limit is a fixed scale, µ jet or µ soft , chosen to best approximate the combination of the hard-scattering functions with parton distributions, although other choices are possible; see Ref. [50]. The conventional choice for many direct QCD analyses has been, for example, µ s = √ŝ /N for the in and fragmentation jets, which resums all logarithms ofN in moment space, to any fixed order in the coupling. The integral now captures all logarithms ofN , but has a branch cut beginning atN = ( √ŝ /Λ QCD ) in the complex N plane. Nevertheless, we can use a minimal procedure to invert the transform [22], by evaluating a contour in moment space that is to the left of the branch cut. This procedure keeps all the logs, at the cost of unphysical, but power suppressed, contributions [23]. The anomalous dimensions, of course, do not depend on this solution to the evolution equation for our jet and soft functions. We give the application of this method to incoming jets first, closely following the discussion of Ref. [21], and then turn to recoil and fragmentation jets, which require a different treatment. For definiteness, we provide solutions that organize all lnN dependence, and turn later to a discussion of other choices of hard and soft scales. in in Eq. (24) are normalized in terms of threshold resummation for the Drell-Yan cross section, which involves only initial-state jets at partonic threshold. When taken at a common Mellin moment N , the squares of theJ (i) in , i = a, b give the threshold-resummed Drell-Yan and Higgs production hard-scattering functions to leading power in N . We define them as a function of a single scale, again labeled µ rf here, as in terms of W (i) DY , the vacuum expectation value of the "cusp" color singlet product of Wilson lines in the color representation of parton i. For incoming jets, the moment variable N i is defined in Eq. (25). Equating the factorization and refactorization scales, µ F = µ rf , this function obeys the renormalization group equation Eq. (27). Its evolution has been discussed widely [43,45,46,18].
To be explicit, the single-log anomalous dimension in Eq. (27) is given by [43] γ J (i) (α s ) = 1 2 Because of (38), this is one-half of the full anomalous dimension associated with the function W We will return to the use of fixed scales below, but here we develop the resummed expressions of direct QCD, in which the infrared scale is chosen to generate all logarithmic dependence in the moment variable, N . We carry out the resummation for µ F = µ rf , choosing µ 2 = √ŝ /N i as the infrared scale and leaving µ 1 = µ rf ∼ √ŝ as the ultraviolet scale in (40). The solution of Eq. (40) then becomes We then use Eq. (28) to reintroduce an independent factorization-scale dependence by extending the integration limit for the integral over the "cusp" anomalous dimension A i , In this form, the prefactor still generates logarithms of N through the running coupling, but we can promote this dependence to the exponent [45], thereby generating all NNLL logarithms from the resulting anomalous dimension. First, we introduce the notation, where the explicit value given in the second line can be found in [21,43], for example. For the expansions of A i and all other functions, we will use the conventions We now note that This enables us to write the resummed in-jet, Eq. (42) in a form where all N dependence is generated in the exponent, Compared to (42), the coupling in the prefactor is evaluated at scale µ rf ∼ √ŝ , and we introduce a new function,D i defined bŷ in terms of the QCD β-function, Eq. (35). The functionD i , which differs from γ J (i) in Eq. (39) by the term proportional to β 0 , also begins at order at α 2 s , Interestingly, the ζ(2) term in γ J (i) is fully cancelled by exponentiating the N -dependence of the jet prefactor. Appendix A collects the explicit low-order expansions of all other anomalous dimensions needed for the NNLL resummed jet functionsJ For completeness, we go on to link the solution, Eq. (46) to another standard expression, with an explicit Mellin transform in the exponent. This form is particularly natural when there is only a single hard scale, as in threshold resummation for electroweak annihilation, or for singleparticle inclusive cross sections in electron-positron annihilation. Details of the transformation to all logarithmic accuracy were given in [51] and reviewed in [45] and elsewhere. For this discussion, we identify µ rf = √ŝ . We then define two shifted R and D functions, and In these terms, and staying at the level of NNLL, we havẽ This expression, with a moment integral in the exponent, is discussed in detail in Ref. [51] for Higgs production (see particularly appendix A there); the only difference here is in the factor R i (α s ), which we normalize to Drell-Yan cross sections.

Recoil and fragmentation jets
The recoil jet and fragmentation jet functions in the factorized expressions Eq. (20) or (24) can be extracted from the singular z → 1 behavior of single-inclusive cross sections in electroweak annihilation. We will describe this procedure, giving some details of how the refactorization that leads to (20) and (24) can be carried out.
The cross sections in electroweak annihilation, e + e − → hX, can be put in factorized form starting from Eq. (5), simply replacing the parton distributions by two delta functions, δ(1 − x a,b ), giving, For ease of comparison to the general form in Eq. (20), we have kept the rapidity dependence relative to the beam direction. For e + e − annihilation, z,η = η and x T are related by z = x T coshη and z c = x T /x T . Equation (52) applies to cross sections for both hadrons and partons, and to analyze the partonic hard scattering we consider the cross section for an observed parton c in the (dimensionally) infrared-regulated theory.
After subtraction of collinear poles from radiated and virtual partons in the fragmentation direction, the inclusive hard-scattering function ω e + e − →c in (52) is infrared finite for z c > z, but with singularities in the limit z/z c → 1. It can be a function only of the remaining invariant mass of the radiated state, which is given byŝ(1 − z/z c ) (note that for e + e − collisionsŝ = s). The hard-scattering function for this cross section then (re-)factorizes at partonic threshold into a short-distance function and two jets, one in the direction of the observed parton and the other of its recoiling partner, with additional wide-angle radiation. As we shall see, for this particular case we will not need a soft function [52].
As an intermediate step in bringing out this structure, we refactorize the inclusive hardscattering function ω e + e − →c into a true short-distance non-radiative factor, denoted H e + e − →c and the function that contains all soft and collinear radiation, denoted Σc, (7)). This is analogous to Eq. (20), but without separation of jets in the final state. As usual, the function Σc will itself require renormalization, so that an additional "refactorization scale", labelled µ rf as above will appear in both it and in the short distance "hard part", H e + e − →c . All y-dependence, associated with the final states, is in the function Σc, and at leading power in 1/ȳ, the function H e + e − →c has contributions only from virtual corrections. Theȳ dependence in Σc links its evolution to that of the fragmentation function D h c (z c , µ F ) in Eq. (52). At the same time, H e + e − →c has a residual dependence on the factorization scale µ F , which will be determined as in the general case, Eq. (36). We are going to give a further factorization of the function Σc to define normalizations for the recoil and fragmentation jets introduced above.
The virtual hard function H e + e − →c in (53) is normalized so that at lowest order the final-state function Σc is a simple delta function setting y to one, At higher orders, the function Σc contains all long-distance dynamics, associated with soft gluons and collinear enhancements in the recoil and fragmentation directions. In the threshold limit, all radiation is forced to be soft, except for particles emitted in the recoil jet direction. In electroweak annihilation, the flavor associated with this jet is naturallyc, as indicated by the subscript on Σ. The singular, long-distance behavior in ω e + e − →c near partonic threshold is given in terms of a sum over states consisting of all radiation except for the observed parton, as produced by the field of the other high energy parton, which initiates the recoil jet. This field, which we will denote bȳ φ r , is joined locally to a lightlike Wilson line that extends to infinity in the direction opposite to the recoiling jet.
To construct the function Σc in a manner that reproduces all singular behavior, we sum over a phase space that matches the phase space of the 1PI cross section near partonic threshold, including all singular regions, and is weighted by squared amplitudes that generate all soft and collinear enhancements. For this analysis, we take equal scales, µ F = µ rf ∼ √ŝ . As for the initialstate jets, independent factorization scale dependence will be reintroduced in the cross section by standard evolution equations. For arbitrary flavor,c ≡ r, the function Σ r is now given in terms of a single scale by [53,54] Σ r ȳ,ȳŝ with the normalization factor N 0 chosen to set Σ r to δ cr δ(ȳ) at zeroth order. In (55),φ r is the partonic field that produces the recoil jet; the trace is over color indices and spin. The familiar As defined in (55), the renormalization of Σ r includes the removal of virtual collinear singularities associated with the presence of the Wilson line, in addition to the collinear singularties resulting from real radiation. This renormalization is equivalent to the subtraction of counterterms defined by the perturbative fragmentation function in (52).
The function Σ r in (55) is the imaginary part of a renormalized "jet function", as defined in Ref. [53] with a light-like Wilson line (and in Ref. [54] with a Wilson line off the light cone). Because only collinear poles have been subtracted in the factorized expression (52), the renormalized function Σ r can still have singular plus distributions associated with radiation collinear to ther direction. For 1PI cross sections in e + e − annnihilation this is the direction of the observed particle (r = c at partonic threshold). This will motivate us to make a further factorization below, into recoil and fragmentation jet functions, as in the general form of Eq. (20).
The sum over states in Eq. (55) matches the phase space for radiation in the limit of y → 1 for single-particle annihilation with a parton observed in directionβ r . In the center of mass frame defined by the timelike momentum p ξ and a lightlike vector p μ r = ŝ/2β µ r , states |ξ consist of soft radiation in theβ r direction and collinear radiation in the direction β r . Note thatŝ enters only in the combinationȳŝ. Collinear-soft radiation is generated by the path-ordered exponential in representationr with constant velocityβ r , originating at the origin, as above. In Eq. (55) applied to e + e − 1PI cross sections, Φβ r plays the dual roles of representing the outgoing observed parton and of serving as a "gauge link" for the partonic field, φ r , rendering the matrix elements in Eq. (55) gauge invariant.
At this stage, it is convenient to take moments, under which the 1PI convolution in Eq. (52) breaks into a simple product of the moments of the parton-to-parton fragmentation function (taking h = c in Eq. (52)) times the moments of Σ r , At leading power, Σ r has an overall behavior of 1/ȳ =ŝ/p 2 ξ times logarithms ofȳŝ/µ 2 rf , so that, as usual in threshold resummation, the function has one less argument in moment space.
We will eventually refactorize the functionΣ r into perturbative recoil and fragmention jet functions, as in Eq. (20). This procedure is based on the evolution equation satisfied byΣ r , which is known. BecauseΣ r is built from the composite parton-Wilson line vertex in Eq. (55), with µ F = µ rf , it again obeys an evolution equation [55,53] of the general "cusp" form (27), with another single-logarithmic anomalous dimension, which we label γ The value of γ (r) Σ reflects renormalization of the composite operator linking the parton field φ r with Wilson line Φ (r) βr in (55). We will determine it by comparison precisely to the threshold-resummed single-particle annihilation cross section, Eq. (52) [52].
The solution to Eq. (58) organizes logarithms of N as an exponentiated integral from scalê s/N to a hard scale, µ rf ∼ √ŝ , We now go through the same steps as for the incoming jets of the previous subsection, starting by promoting the N -dependence of the coupling of the evolved coefficient function to the exponent where we introduce factorization scale dependence, as in Eq. (42) for the in jets, and where we defineB by analogy to the shift in the functionD for initial state jets (see Eq. (47)).
At this point, we can determineB and then γ Σ from the literature. To do so, we further rewrite (59) with moments explicit in the exponent, again as for incoming jets, Reexpressing the exponent in the explicit Mellin form leads, as above, to a new prefactor at the hard scale, as shown explicitly, and to a shift fromB r to B r . The latter depends on the kinematic range µ 2 rf < µ 2 < (1 − y)ŝ in (62), and turns out to be one-eighth of the corresponding shift in the "Drell-Yan" function D(α s ), which appeared in the discussion of incoming jets above. Again following the steps taken in Refs. [51] and [45], we find By taking µ 2 rf =ŝ, we can compare to the explicit expression for B r as it appears in single-particle inclusive resummation at center-of-mass energyŝ, as given [52]. At O(α 2 s ), we find forB r , This is the anomalous dimension we will use below. We recall the values of the B (2) r in (63) as given in the literature [52] in Appendix A.
We can now solve Eq. (61) to find an explicit expression for γ (r) Σ (α s ). This requires the one-loop Mellin space prefactorsΣ r (1, 1, α s ), which we can compute directly (see Sec. 4.3 below) to find Using these explicit forms in Eq. (61), we have for the two-loop term of γ At this point we may recall the relation observed in Ref. [56,57], valid to two loops, where P q,δ is the coefficient of δ(1−x) in the quark DGLAP splitting function (see App. A), and as above D i (see Eq. (50)) is the single-logarithmic anomalous dimension for Drell-Yan resummation [3,4,56,51]. In terms of the anomalous dimensions γ J (q) (α s ) of Eq. (39) and the above γ q Σ (α s ), relation (67) becomes A similar result holds for the gluon. We also note that at order α s , the single-loop anomalous dimension is given by the delta function piece of the splitting function, which is consistent with jet anomalous dimensions in direct QCD and soft-collinear effective theory.
In the form of Eq. (59), the integral over scales µ that appear in the running coupling can be reorganized in a trivial fashion to separate the recoil jet from the fragmentation jet, as presented, for example, in Refs. [56,57], where as shown, we introduce factorization scale dependence into the fragmentation jet as we did for the incoming jets. Note that dependence on the ratioŝ/N 2 µ 2 rf cancels in the product of the recoil and fragmentation jets.
The factorization in Eq. (69) is, of course, not unique. For application to QCD scattering, we normalize the outgoing jets so that they preserve the one-loop structure of the soft anomalous dimension matrix at two loops [58]. We can do this by choosing the fragmentation jet to be equal to the incoming jets of Drell-Yan as a function of N . This is the natural choice for doubleinclusive annihilation, as noted in [59]. It is also the manner in which the soft anomalous dimension matrices were normalized and computed in [58]. The formal definition of the eikonal fragmentation jet functions differs from that for the incoming jet function in being defined by the eikonal doubleinclusive annihilation cross section, rather than the eikonal Drell-Yan cross section. Soft radiation is thus emitted from an outgoing color-singlet pair of Wilson lines, rather than the incoming color singlet pair ofJ (i) in , Eq. (38). The equality of such eikonal fragmentation functions and parton distributions defined at fixed energy is shown in Ref. [59].
With this normalization, the fragmentation jet is given by the expression for incoming jets, Eq. (46), involving both the cusp anomalous dimension and the function D. We will choose partonic subscript c for the fragmentation jet to emphasize that our choice here extends beyond e + e − annihilation, This definition is the same as in Ref. [21].
Using Eq. (62) forΣ r and (70) for the fragmentation jet in (69), we find that the recoil jet includes the fullB r term in the exponent, while the terms involving A r (α s ) change, We see here thatJ (r) rec obeys the evolution equation given in (29), with a single-log anomalous dimension that is a combination of the functionsB r andD r . We note that although the function D r cancels in the 1PI electroweak annihilation cross section, these expressions can be chosen as the factorizing jets in QCD cross sections, where the outgoing partons are not necessarily of the same flavor. Alternative choices for the jets are possible, but in general will lead to soft anomalous dimensions with additional terms proportional to the identity matrix at two loops. Although the separation of recoil from fragmentation in Eqs. (70) and (71) is a simple reshuffling of the integral, it has a direct physical interpretation because of the relation of the scale that appears in A r (α s ) to the momenta of "web functions", which, as discussed for example in Ref. [60], reflect the transverse momentum and hence direction of soft radiation. Closed expressions for NNLL expansions of the jet functions in Eqs. (70) and (71) are again given in Appendix B.
In summary, for the short-distance coefficient function for electroweak annihilation, we reexpress Eq. (53) as a special case of the general factorized cross section, which exhibits the separation of the fragmenting and recoil jets in the simplest case, where the µ F dependence of the function H is determined by Eq. (36), in this case with only a single term, for the fragmenting parton. It is straightforward to verify the consistency of this expression with the explicit results for single-inclusive annihilation given in Ref. [52].

One-loop partonic jets
For later reference, we now expand each of the jet functions to order α s and transform the result back to z space. For each, we write the result in the form For the initial-state functions discussed in Sec. 4.1, we use the explicit moments given by Eq. (25).
If the underlying Born process is ab → cr, with parton c fragmenting, we have where In Eq. (74) we have neglected contributions that are suppressed near threshold. Note that we have also suppressed the dependence on the factorization scale, which may be easily reconstructed.
The construction of jet functions for the outgoing fragmenting partons was presented in the previous subsection. In summary, we compute the renormalized function Σ r for each parton flavor, and separate the fragmentation and recoil jets as in Eqs. (62)- (70). The Σ r for quarks and gluons are not related to each other by simple exchange C F ↔ C A . At one loop, we find for quarks and gluons from (55) The jet function associated with the fragmenting parton is given by Eq. (70), and is found at one loop by a simple expansion and inverse transform, The jet functions for the recoiling quarks and gluons are defined as the remaining part of the final-state function Σ r , Eq. (55), and are found from (76) by removal of the fragmenting jet (77), and Here we have introduced b 0 = β 0 /(4π) = (11C A − 2N f )/(12π). We note that the logarithmic terms are the same as those found from the expansion of the resummed expression for the recoil jets, Eq. (71).
It is useful here to make again contact to e + e − → hX, which we used in Sec. 4.2 to define the outgoing jets. The hard function in this case is obtained from the known [61] one-loop virtual correction: where we have suppressed the label e + e − → qq for the process, and where We thus find where ω e + e − →c,(1) is the one-loop inclusive hard-scattering function of Eq. (52), calculated to NLO in [62]. As designed, our formalism thus reproduces the full NLO correction near partonic threshold, without the need for any additional soft function. For general single-particle cross sections in hadronic scattering, a soft function is of course necessary.

Soft Functions for Single-Particle Cross Sections
Soft radiation in the central region between the jets is well approximated by light-like Wilson lines in the color representations of the partons participating in the hard scattering [39,26,22,19]. We begin our discussion of the necessary soft functions by briefly reviewing their resummation.

Resummation for soft functions
The soft functions in the refactorized expressions, Eqs. (20) and (24) are matrices in the space of color exchange tensors for the partonic process ab → cr [39,7]. We will specify their definitions in the next subsection, and here introduce some notation, since we will use both momentum and moment space. The soft function scales as an overall w −1 S , with additional logarithms in the combination w S √ŝ /µ rf . The soft function is constructed to have at most one such logarithm per loop. Higher logarithms are associated with collinear enhancements, which are universal, and are factored into the jet functions [26,7]. To identify fixed-order terms, we introduce the momentum-space notation, where here and in the following we usually suppress the explicit labeling of the underlying partonic process, ab → cr. Note that from the refactorization expression, Eq. (20) the argument of each S (n) i appears in a convolution, whose range is 0 ≤ w S ≤ŝ 4 /ŝ. The corresponding moment-space expansion of the last equation will be denoted bỹ where as beforeN ≡ N e γ E . In the following section, we will provide field-theoretic definitions of these soft functions, and give an example of the calculation ofS itself provides a series of NNLL logarithms, as we now review, in addition toS (1) 1 , which is proportional to the anomalous dimension matrix Γ, and contributes at NLL.
The resummed soft factor in moment space is the solution to Eq. (30) [39,26,22] in terms of the refactorization scale, (85) The second factor on the right-hand side of this expression is the soft function in moment space, Eq. (84), with a scale choice that makes all its logarithmic terms vanish. A full NNLL resummation takes into account logarithms due to the expansion of the running coupling for the one-loop soft function at scale √ŝ /N , in much the same way as it includes logarithms from the function R i (1, α s ) from incoming and fragmentation jets, as in Eqs. (42), (45) and fromΣ r (1, 1, α s ) from the recoil jet, in Eq. (71). This factor serves as the boundary condition for the evolution of the soft function from infrared to ultraviolet scales. The factor that resums the evolution logarithms of the moment variable in Eq. (85) is given by the ordered exponential [7] S ab→cr N , α s (µ rf ), with P denoting path ordering. The process-dependent anomalous dimension matrix, Γ ab→cr is determined entirely from virtual corrections. As observed in Ref. [58], the two-loop anomalous dimension matrix, Γ ab→cr is proportional to the one-loop matrix, with K = C A (67/18 − π 2 /6) − 5N f /9. This specific relation does not extend to three loops, as was demonstrated by explicit calculation in Ref. [63], but this does not lead to qualitative differences in the analysis beyond NNLL. As in Refs. [20,21], exponentiation of these matrices is readily carried out numerically, by iterating the exponential series to an adequately high order.

Operator definitions
As in the dihadron case, the elementsS LI of the soft matrix of Eq. (24) are computed using the method described in Ref. [26], which, however, we must adapt in a significant way to single-particle inclusive cross sections. The all-orders form is clearest in moment space, where it is given as the ratio of the moments of a fully eikonal cross sectionσ ab→cr LI and four factorized jets, two to absorb the factorizing collinear singularities of the incoming parton lines, and two to absorb the collinear singularities of outgoing lines, all in eikonal approximation: where as aboveη is the rapidity of the fragmenting jet in the partonic center of mass frame, and s sets the scale of the invariant mass of the partonic system in the single-particle inclusive case we are considering. As indicated, dimensional regularization with D = 4 − 2ε dimensions is used to regulate divergences.
In order to formulate operator definitions for the quantities on the right-hand side of Eq. (88), we introduce a slight generalization of the definition of the Wilson line in Eq. (56), which will appear both for the numerator and for the "eikonal jet" functions in the denominator of (88). For 2 → 2 scattering, the ends of two incoming and two outgoing Wilson lines are coupled locally by a constant color tensor C I , and we define These operators will produce radiation at all scales and directions, including collinear radiation. As described in Refs. [39,26] and below, the incoming and fragmentation jets,j in andj fr , respectively are constructed to match collinear singularities and radiation phase space in the partonic threshold limit, avoiding double counting with the partonic jet functions in the re-factorized hard function, Eqs. (20) and (24). In the same way, we also define singlet operators that link two lines in conjugate representations, extending from the infinite past and joined at the origin by a color singlet tensor: We will use these to construct the incoming eikonal jets.
In these terms, for the incoming eikonal jets, we construct the eikonal analogs of partonic (Drell-Yan) annihilation. Unlike the case of the final state jets below, the phase space for the initial state jets is defined by a total energy, and is hence finite. The kinematics of the process is reflected in the rescaled Mellin moment variables, as in Eq. (25) [7]. The "in" jets constructed in this way are found in Ref. [21] on dihadron production, and are given by With this choice, j (a) in 2 is the eikonal Drell-Yan cross section, computed at two loops in [43]. In fact, the eikonal jets that remove collinear singularities from the eikonal cross section are the same as the incoming jet functions that appear in the re-factorized hard-scattering functions, defined as in Eq. (38). This is simply a reflection of the multiplicative nature of factorization in moment space.
For the outgoing jets, we turn to eikonal single-particle inclusive e + e − annihilation. Here, the eikonal cross section is defined at fixed values of the invariant mass of all radiation recoiling against the observed particle [7,8,17], as in the definition of the recoil and fragmentation jets, derived above from functions Σ r , Eq. (55). In particular, because our cross section is defined at fixed momentum fraction y, a light cone fraction in the direction of the observed momentum p c , we must incorporate the limitation on the energy of radiation in the direction of the recoil jet in the definition of the eikonal cross section and its eikonal jet subtractions. This is because fixing y alone does not result in a finite phase space inσ ab→cr LI . Specifically, fixing y and hence s 4 still allows collinear divergences in the recoil direction from arbitrarily large energies. The situation is to be contrasted to pair invariant mass threshold resummation, where a fixed energy automatically imposes a limited phase space. In the present case, we must truncate the sum over radiation collinear to the recoil direction.
To this end, we match the physical phase space for partons in the soft function to the phase space near 1PI partonic threshold to cut off unbounded collinear behavior. The details of the collinear truncation will cancel in the ratio of Eq. (88), because collinear partons factor from soft gluons emitted at fixed angles both in the numerator and denominator. In defining the space of states over which to sum, we can thus replace the partonic recoiling jet by a single particle, whose momentum we will denote by p where we integrate over the n-particle phase space, requiring p (h) 4 2 = 0.
In the p c /p R c.m. frame, with p c in the plus direction, p R has only plus and minus components, which we denote as where in the second equality β c is the lightlike vector in the direction of the observed particle, and β r is the "opposite-moving" lightlike vector in the direction of the recoiling jet at partonic threshold in this frame, with β c · β r = 1. This definition will enable us to evaluate integrals over partonic phase space in other frames, and this is the form we will use below.
In terms of the operators w I , the eikonal cross section is defined by sums over states |ξ = |{k i } radiated freely by the Wilson lines, subject only to the momentum constraint of Eq. (93) involving p R , Eq. (94): where p ξ is the momentum of state |ξ . Thus we define the eikonal soft function with exact eikonal matrix elements, but as a sum over the states |ξ , consisting of radiated gluons, and at high orders, quark pairs. At NLO, for the soft function, S (1) , we need only the single-gluon final state, n = 1, for which the condition (93) is equivalent to with k the momentum of the radiated gluon. Again, the momentum p is introduced only to define the phase space; in the evaluation of the eikonal cross section at one loop below, only the recoilless "eikonal" vector p 4 and the total recoil momentum p R appear.
The definition of the eikonal fragmentation jets is designed to match the phase space of the partonic recoil jets. For β c the fragmentation direction, it is given by the same eikonal fragmentation function as for the full jet functions in Eq. (70) [21] Similarly, the eikonal recoil jet is the eikonal analog of the partonic recoil jet extracted from Eq. (55), using the same phase space as in the eikonal cross section (95) and is given bỹ fr the same MS distribution as in Eq. (97), whose flavor is defined by the parton initiating the recoil jet.

One-loop soft functions
We are now ready to determine the finite soft function at one loop and beyond, taking the simplest case of the scattering of quarks of different flavors. The calculation of soft matrices for the other partonic processes follows exactly the same pattern. As we discussed above, and as we shall see explicitly, the soft matrices for single-particle inclusive cross sections reflect the phase space of this process and differ from related cross sections with the same partonic channels, such as dihadron production, treated in a similar fashion in [21].
The explicit calculation of (S ab→cr ) LI at one loop as given here is equivalent to the procedure described in Ref. [21]. The functions on the right side of (88) are are written to one loop aŝ is the moment-space tree-level soft matrix, Eq. (84) corresponding to the particular partonic process [39], and wherej can be any of our eikonal jet functionsj in ,j frag ,j rec . In moment space, the collinear singularities of the eikonal cross sectionσ LI are cancelled by those of the incoming and outgoing jet functions, constructed as above. Expanding the soft function to first order, we have in moment space The eikonal cross section at order α s may be determined from the phase space integrals (defined here in Feynman gauge) with the momentum p R defined in Eq. (94),ȳ ≡ 1 − y, and the kinematic variable v as defined in Eq. (75). These integrals appear times color coefficients, which we label as R ij , that are color tensors and depend on which Wilson lines, i, j = 1, . . . , 4 are connected by the emitted gluon.
Similar integrals are encountered in the treatment of threshold resummation for dihadron cross sections [21], which extend over a different phase space.
The momentum-space eikonal integrals I ij are, of course, the same in any representation for the Wilson lines, and are given in D = 4 − 2ε dimensions by We note here that these integrals involve only logarithms of the kinematic variable v, rather than the dilogarithms found for dihadron threshold resummation [21]. The explicit forms found here are, of course, necessary to reproduce singular threshold behavior at one loop and beyond.
We now illustrate the calculation of the soft matrix using these ingredients, specializing to our example for the partonic channel with quarks of different flavor, qq → qq . In this case, the Wilson lines are in the fundamental representation, and the color-space matrices R ij describe how single-gluon exchange mixes the couplings of the Wilson lines that represent soft radiation [26]. For definiteness, we choose a basis of t-channel color exchange between quarks of distinct flavors. The lowest-order soft matrix is independent ofη and in this basis given by [27] We note that the S With lines in the fundamental representation, the R ij mix singlet and octet exchange. Acting on the amplitude represented by vectors with the t-channel color singlet coupling of the Wilson lines in the first entry and the t-channel octet in the second, they are given by In these eikonal factors, the interference between initial-and final-state emission has a relative minus sign which we exhibit here, changing the notation slightly from that of Ref. [21]. Notice that the sum of the R ij is proportional to S (0) . A corresponding result holds for all channels, and follows from the gauge theory Ward identities. This ensures that at one loop double poles factorize and cancel in Eq. (100).
Together with the dI ij /dȳ, these matrices define the eikonal cross section (95) at one loop. Taking moments, we find where as before N a = N v, N b = N (1 − v). On the right-hand side of this relation, we have suppressed the subscript denoting the process, qq → qq for the soft functionS (0) 0 and the anomalous dimension matrix Γ (1) , which is the first-order term of the matrix defined by Eq. (86). For the process at hand (and in fact for all quark-quark scattering processes), we have [27] Γ qq →qq , (1) where T = ln(1 − v) + iπ, U = ln(v) + iπ.
Following Eq. (88), to derive the soft function, we next need to divideσ qq →qq LI by the eikonal jet functions. In computing the ratio, we need consider only real-gluon contributions. Virtual corrections to these scaleless integrals can be defined as pure counterterms, which cancel the infrared poles of the real contributions, as in Ref. [64] for example. All double poles also cancel in the ratio, leaving only a single pole, which can be considered part of the renormalization of the soft function. To order α s the real-gluon contributions to the incoming eikonal jets defined in (92) are, in momentum (y) space, with the kinematic variable v defined as in Eq. (75), and as usual, For the observed (fragmenting) parton, defined as in Eq. (97), we have at order α s for the real-gluon contributions with c = q or g. Finally, to obtain the one-loop jet function for the unobserved recoiling parton, r as defined in (98), we subtractj with r = q or g. Taking moments and dividing out the jet functions from (105) cancels all collinear and infrared poles, and we obtain at O(α s ) After MS renormalization, the terms proportional to the anomalous dimension matrix result from the evolution of the zeroth order soft function, and the remaining expression is the one-loop finite term of the soft function, as in Eq. (84): for qq → qq . According to Eq. (85), to NNLL accuracy and beyond, this function will appear in the resummed hard scattering, Eq. (24) times α s evaluated at scale √ŝ /N . We note that the first-order soft matrix may be cast into the form where C a = C b = C c = C r = C F , and with R 12 as in Eq. (104).
The soft matrix constructed here, along with its matrix anomalous dimension, bears a close relation to the "wide-angle" soft function computed in the context of top production in Ref. [19]. In fact, non-diagonal entries in the soft function are the same (in the same basis), but the soft matrices are not identical. Our evolved soft function acquires only a single logarithm of the moment variable per loop, while the soft function in [19] has up to two logarithms per loop. These double logs are color-transfer independent, however, and the difference is due to the different choices in soft-collinear factorization commonly made in SCET, compared to direct QCD. This comparison is discussed, for example, in Ref. [47]. In particular, in the direct QCD treatment chosen here, all double logarithmic behavior is factorized into jet functions. The consistency of the non-trivial single-logarithmic evolution is another confirmation of the underlying consistency of the two treatments.
Soft functions for all partonic subprocesses are constructed in the same fashion, starting from eikonal cross sections, and dividing out eikonal jets. Results are presented in Appendix C. We found that in each case, the result takes the form given by Eq. (112), except that for processes with a qq ( ) initial or final state, the sign of the last term needs to be reversed. 6 The Resummed Short-Distance Function and Moment Inversion

The resummed inclusive hard-scattering function in moment space
We are now ready to combine our previous results and present the resummed hard-scattering function, Eq. (24) in transform space, to leading power in the transform variable N .
We will base inverse transforms of the expressions above on a hard scattering function written almost entirely in terms of logarithmic integrals. Indeed, this is the form on which our derivation in Sec. 4 of the resummation is based. The expressions we need are Eqs. (43) and (46)  We give the result for a specific choice of refactorization scale, µ rf = √ŝ . This is a scale that simplifies existing expressions for the hard-scattering function [42]. We comment below on alternative refactorization scale choices.
In these terms, the full expression with µ rf = √ŝ , to NNLL, is a product of exponentials associated with the jets, multiplied by a trace that links the hard and resummed soft matrices, Here, the kinematically-rescaled moment variables N a and N b for the incoming partons are given by N a = N (−û/ŝ), N b = N (−t/ŝ) (see Eq. (25)), and we define N c ≡ N .
In practical applications of the result in Eq. (113) we will typically need an expansion in terms of the coupling at a fixed scale, µ R . For instance, matching to fixed-order calculations, we will write the hard function in the form in which the refactorization scale appears in coefficients. These coefficients have up to two logs in µ rf / √ŝ per loop, because of the renormalization group equation satisfied by the hard function, Eq. (34). Additional single-logarithmic terms appear from the reexpansion of the running coupling. Dependence on the renormalization scale µ R now begins, as usual, at the next uncalculated order. As noted above, the factorization scale dependence of the hard matrix is known order-by-order through Eq. (37). For our choice of µ rf = √ŝ , we may use the notation, The general form of these hard functions is, of course, similar to that given in Ref. [21] for dihadron cross sections, but with different final-state jets and soft functions, as developed above. Specialized to the case of prompt photon production (c = γ), these results extend the NLL treatments given previously in the literature [8,10,11]. At NNLL (and beyond) they are also consistent with the soft-collinear analysis for prompt photons in Ref. [18], as we discuss in Sec. 6.3 below, with similarities in moment space as explored for Drell-Yan in Refs. [65,45].
In much the same way as Eq. (114), the resummed exponents in Eq. (113) will also start to depend on a renormalization scale µ R in fixed-order expansions. In fact, even when keeping the all-order resummed form of Eq. (113), dependence on a scale µ R will arise in the usual "minimal expansions" [22] of the exponents, despite the fact that the full exponents are strictly independent of such a scale. This is because the perturbative expansion of the resummed exponents to a desired logarithmic accuracy necessarily truncates the perturbative series. The dependence of the resummed exponents on µ R is explicitly seen in the expansions given in Eqs. (142) and (147) in App. B. It is also the dependence that we will explore numerically in Sec. 6.5.
Equation (113) is the central result of this paper. As observed above, we have given it for a specific choice of refactorization scale, µ rf = √ŝ , because this choice simplifies logarithms in the hard function, H. Other choices are, of course, possible, leading to changes in the expression, and presumably, in numerical results. The simplest modification is to vary µ rf by a factor around √ŝ , say µ rf = ζ √ŝ , with 1/2 ≤ ζ ≤ 2. In Eq. (113), this can be implemented by replacing √ŝ by ζ √ŝ wherever it appears explicitly without a factor ofN . Logarithms of ζ will also appear in the explicit expansion of the hard matrix H. Of course, variations can also be implemented independently in both the hard and soft endpoints of each scale integration in Eq. (113) [50].  v). For the purposes of this resummation formalism, none of these terms is taken to be parametrically large, although this depends to some extent on the process-and energy-dependent kinematics. A formalism with applicability to large partonic rapidity, as treated in electroweak annihilation in Ref. [66], will require further development for QCD hard scattering.

Comparison to NLO
With the complete short-distance functions in hand, we have all the ingredients necessary to compare to the full one-loop calculations available in the literature [5,6]. Expanding the resummed partonic cross sections to NLO, we should recover all leading contributions that arise near partonic threshold in the full NLO cross sections. These are all terms with distributions of the form (ln(ȳ)/ȳ) + , 1/ȳ + , and δ(ȳ) where, as before,ȳ = 1 − y. This is a powerful test of our resummation procedure. It is worth pointing out that we do not actually use the NLO cross section to determine the δ(ȳ) pieces in our resummed cross section, but that we independently predict these pieces from our expressions given above. We note that in [21] we showed that the resummation formalism fully reproduces also the NLO hadron pair cross section near threshold. We now extend the comparison to single-inclusive production. This primarily tests our final-state recoil jet function and our new one-loop soft function for this case. We use the simplest process qq → qq as an example and present the corresponding results for all other partonic channels in Appendix D.
Expanding the inclusive hard-scattering function in Eq. (24) to first order we have, again transforming to momentum space, and using the form of the one-loop soft function, Eq.
where on the right we have dropped the label ab → cr of the subprocess. We now present explicit results for the various contributions in the above equation for the process qq → qq . Using the results in (74)- (79) and in (103), (111), and defining we have: and Tr With these, Eq. (116) reproduces the full NLO [5,6] for qq → qq near threshold, including all δ(ȳ) contributions. Appendix D collects the expansions for all other partonic processes; we have checked that in each case the resummed formulas reproduce NLO near threshold.

Resummed cross section with N -independent scales
As noted above, in Ref. [18] and other effective theory resummation studies, it is customary to use fixed, rather than N -dependent, infrared and ultraviolet scales in solutions to the evolution equations that generate threshold resummation. We can do this here as well, which leads to an expression that generalizes the results of Ref. [18] for prompt photons to final-state hadrons. Following this approach, we replace √ŝ /N by a soft scale, µ s for the incoming and fragmenting jets and the soft function, and ŝ/N by a jet scale, µ j for the recoil jet. In the solutions, dependence on the moment variable N is retained in the prefactors for the various jet functions, so that the jet anomalous dimensions γ J and γ Σ appear in the exponents rather than the alternative anomalous dimensions D and B.
The solutions for the incoming and fragmentation jets can be found from Eq. (40), and then the recoil jet from (59) and (65). Using these for the resummed cross section we obtain, instead of Eq. (113), For N -independent choices of µ j and µ s , logarithms of N do not all appear from the exponents [45]. Again, we recall thatN c =N here.
For compactness and ease of comparison to the literature, we introduce notation in the spirit, if not the exact letter, of that found in the treatment of prompt photon cross sections in Ref. [18], where the final definition applies to γ = γ (r) Σ , γ J (i) . (As a notational point, the "Sudakov" factor S(µ 1 , µ 2 ) should not be confused with the soft function.) In this notation, the moment-space resummed hard-scattering function Eq. (120) becomes When specialized to prompt photon production, this result is consistent with Ref. [18], with a different organization of logarithms in the moment variable compared to Eq. (113). We now turn to inversions of these transforms to find the leading threshold behavior.

The Mellin inverse
To apply any of the moment-space resummations given above for phenomenological applications, we must perform a Mellin inverse. The most direct approach is found from Eq. (12), with a suitable contour C. This result is then convoluted with the parton distributions to give (see Eq. (13)) In principle, because the moments of the fragmentation functions fall off quite rapidly, the numerical inversion of the moments will allow the integration against the parton distributions to be carried out numerically. This is analogous to the procedure in Ref. [21] for the threshold resummation of cross sections at fixed invariant mass and rapidity difference, although in that case the roles of the fragmentation functions and parton densities in Eqs. (123) and (124) are reversed. For dihadrons at NNLL, this worked well by defining the inverse transforms like Eq. (123) following the "minimal" definition of Ref. [22].
In the "minimal" procedure, the contour C in the complex N -plane is chosen to cross the real axis between the origin and the branch cut associated with the leftmost Landau singularity, which occurs in the resummed exponent (113) whenN = √ŝ /Λ QCD . The presence of this singularity allows for power-suppressed but nonzero contributions in the unphysical region z > 1, where for dihadrons, z = M 2 H /ŝ, with M H the dihadron mass. These contributions were indeed numerically small for the phenomenologically-relevant kinematics studied in Ref. [21].
Sample evaluations, however, show that single-particle inclusive cross sections are much less stable against contributions with unphysical origin (z > 1), even though they are still powersuppressed. This increased sensitivity can enter because once z is greater than unity, the partonic fractions x a and x b in Eq. (124) can become unphysically small (see Eq. (7)). Through Eq. (2), this leads in turn to very large values of partonic rapidityη, which affects the limits of integration in the resummed exponents, as well as the soft functions, through their anomalous dimensions. We know of no arguments that limit the normalization of such unphysical contributions.
Another approach to using the resummed inclusive hard-scattering function (113) avoids the Landau pole, and unphysical contributions to the inverse transform [67]. In this approach, usually adopted in soft-collinear effective theory treatments of resummation, we replace the N -dependent solutions to the evolution equations for the soft and jet functions by N -independent jet and soft scales: µ j and µ s , as discussed in the previous subsection. Given such choices, the moment inversion ofω ab→c can be done explicitly, using only the identity where C is a again a contour in the N -plane that can now be taken to the right of all singularities, and where in this form there are 1/N or O(1 − z) corrections. In this analysis, ζ represents a sum of integrals over any anomalous dimensions between the chosen soft or jet and the hard scale. Explicit N -dependence in the soft function is accounted for through derivatives with respect to the variable ζ. In the general case, logarithms of N appear in soft and jet functions, and one uses This is readily applied to Eq. (122). Defining the negative of the sum of lnN 2 coefficients in the exponent as and recalling the definitions of theN i in Eq. (25), we find from Eq. (122) In this form, if we set µ F = √ŝ , we can reduce to the case of prompt photon emission, c = γ, and compare to the effective theory treatment of Ref. [18], taking µ = √ŝ in their notation. To do so, we only need the anomalous dimensions associated with the two partonic processes, pair annihilation and QCD Compton scattering, found in Ref. [7]. We find that the two expressions are consistent in explicit dependence on the soft and jet scales, Sudakov factors and other anomalous dimensions, the relevant sums of which are equal at least to two loops.
In closing this discussion, we note that the advantages of both of these well-explored approaches may be incorporated by a hybrid choice, already suggested by the analysis presented in Ref. [68], given, for example, by for the soft scale, with µ larger than Λ QCD . This is just the sum of conventional direct QCD and SCET boundary conditions, differing significantly from a "pure" √ŝ /N choice only when N is very large, where it bottoms out at a "jet" or "soft" scale µ in the SCET language. The "Landau" branch point is now moved to a large negative positionN L = − √ŝ /(µ − Λ QCD ). Qualitatively, a singularity at this position affects the inverse transform to variable z as an additive contribution that scales as z − √ŝ /µ . The influence of this nonperturbative singularity is thus suppressed exponentially. This is more or less the same power structure as the Landau pole in standard "minimal resummation" [22], but now without contributions for z > 1, and with the non-perturbative singularity far from the origin to the left of the contour. We shall not pursue phenomenological implementations of this or other methods to invert the transform of our resummed expressions in momentum space here. We anticipate exploring a formalism using Eq.
(129) and other possibilities in future work. For now we will restrict our numerical analysis to a brief exploration of fixed-order expansions and scale dependence of the resummed cross section, setting aside z > 1 contributions associated with the Landau pole.

Numerical results for fixed-order expansions
We have seen in Sec. 6.2 that the order α s expansion of the resummed cross section reproduces the singular behavior at threshold for each partonic channel. A natural first test is to compare the numerical result of the expanded resummation formula to the full NLO with realistic choices of kinematics and parton distributions and fragmentation functions. In addition, given that resummation provides insight into the size of beyond-NLO effects, we will also explore expansions to NNLO.
In the following, we will consider neutral-pion production pp → π 0 X at two center-of-mass energies, √ s = 31.5 GeV and √ s = 200 GeV. The former corresponds to one of the energies used in the Fermilab fixed-target experiment E706 [69], while the latter is relevant for experiments at RHIC. Although we will not present any actual comparisons to data, we adopt the proper kinematics for these two cases, integrating over pseudorapidity |η| ≤ 0.75 in the E706 case and over |η| ≤ 0.35 for RHIC, corresponding to measurements by the PHENIX collaboration [70], and considering the cross sections as functions of pion transverse momentum, p T . We will use the CT14 parton distribution functions [71] as implemented in the LHAPDF database [72]. CT14 provides both NLO and NNLO sets of parton distributions, which is useful for our expansions. We compute LO and NLO cross sections with the NLO set of parton distributions, and the NNLO expansions with the NNLO ones. For the π 0 fragmentation functions we use the set of [73], which is at NLO. As discussed above, we need Mellin moments of the fragmentation functions, whereas the set of [73] is available as a numerical code in z-space. Technically, we obtain the moments by performing a fit to each of the fragmentation functions for a given set of scales. The functional form of the fit is chosen such that one can easily take its Mellin moments. We have checked that this approach works to about 1% accuracy in the kinematic regimes of interest to us.
As discussed after Eq. (114) and shown by Eqs. (142) and (147) in App. A, the explicit NNLL expansions of the jet functions induce dependence on a renormalization scale µ R . As is usually done, we will choose µ R and µ F of order p T . Figure 1 shows various "K-factors" for pp → π 0 X with E706 kinematics with µ F = µ R = p T , where with k = 1, 2. The crosses show the K-factor corresponding to the full NLO result of Ref. [6]. We see that the NLO K-factor is large, exceeding 2 throughout the p T regime considered. The lower solid line shows the NLO expansion of the resummed cross section. The agreement of the NLO expansion with full NLO is excellent, at about 2% or better. This is an improvement over previous comparisons for the rapidity-integrated cross sections given in Ref. [15]. It provides confidence that resummation indeed captures the dominant contributions to the cross section and motivates the study of higher-order expansions of the resummed cross section as a method of obtaining accurate results beyond NLO. ab→cr , although we include the scale logarithms in the hard function through second order. Evidently, the expansion to higher orders leads to further sizable enhancements, especially at high transverse momentum where the threshold logarithms become more and more sizable. This result is in line with what was found in Ref. [15].
An interesting observation is that, although the threshold logarithms provide a large part of the enhancements seen at NLO, the one-loop hard functions H (1) ab→cr are numerically very important as well. This is shown by the lower dashed line, which again presents the NLO expansion of the resummed cross section, but this time without the contributions by the H (1) ab→cr . Specifically, in the notation of Eq. (115) we set the O(α s ) correction in H ab→cr (α s , 1, 1,η) to zero, but keep the scale logarithms of H ab→cr that arise at that order. Clearly the result is much lower and falls short of the full NLO result. We recall that in Mellin-space the H Given the importance of the δ(ŝ 4 ) contribution for obtaining a good agreement between the NLO expansion of the resummed cross section and full NLO, one may wonder how accurately the NNLO expansion shown by the upper dashed line in Fig. 1 will really match the full NNLO cross section. Among the five "towers" of leading NNLO corrections near threshold, α 2 s ln 3 (ŝ 4 /ŝ)/ŝ 4 + , α 2 s ln 2 (ŝ 4 /ŝ)/ŝ 4 + , α 2 s ln(ŝ 4 /ŝ)/ŝ 4 + , α 2 s 1/ŝ 4 + , and δ(ŝ 4 ), the first four are fully accounted for by our formalism described in the previous sections. The δ(ŝ 4 ) contribution, however, requires knowledge of the two-loop hard functions H (2) ab→cr , as well as of the presently unknown two-loop soft functions S (2) ab→cr and the O(α 2 s ) corrections to the R i andΣ r . Recalling the dominance of the H (1) ab→cr in the NLO δ(ŝ 4 ) contribution seen above, one might expect that the H (2) ab→cr are equally important for the corresponding NNLO terms. Fortunately, the complete set of H (2) ab→cr has been given in Ref. [42], so that we may include it in our studies. The color basis adopted in that paper differs from the one we use, but it is relatively straightforward to transform the results to our basis. Reference [42] also provides the results for the H (1) ab→cr , and we have verified that after the change of basis all our H (1) ab→cr are correctly reproduced.
The upper solid line in Fig. 1 shows the NNLO expansion when the full two-loop terms H ab→cr are for the NLO one, the result shown would include all five leading NNLO terms near threshold and hence be expected to provide a faithful estimate of full NNLO. We observe that the H (2) ab→cr lead to a further significant enhancement of the K-factor. As expected, this enhancement is moderated toward higher p T where the plus distributions become more and more dominant in size. Figure 2 shows similar results, now for pp collisions at RHIC energy √ s = 200 GeV, again using µ F = µ R = p T . Some of the qualitative features from the previous figure carry over: Again the NLO expansion of the resummed cross section matches the full NLO one very accuratelydespite the fact that we are on average further away from partonic threshold due to the higher collision energy. The contribution by the H (1) ab→cr is again crucial in order to achieve this. In fact, the H (1) ab→cr are relatively more important than in the fixed-target case, as one would expect. The NNLO terms in the expansion again provide an enhancement of the cross section; this time the enhancement becomes particularly pronounced only when the H (2) ab→cr are included as well. Figures 3 and 4 show the same calculations as Figs. 1 and 2, respectively, but now choosing scales µ F = µ R = 2p T . We observe that for this scale choice the K-factors turn out to be even larger than the ones we found in the previous figures. The hard functions H (1) ab→cr and H (2) ab→cr turn out to be relatively less dominant, although still important, for this scale choice.
The results we have shown so far suggest a relatively strong scale dependence of the fixed-order expansions. To investigate this further, the upper left part of Fig. 5 shows the µ R dependence of the invariant cross section for E706 kinematics at p T = 5 GeV, keeping a fixed value µ F = 2p T . We vary µ R no further down than p T , in order to avoid having µ F and µ R too different. The lower dashed line shows the LO result, which exhibits a very large scale dependence. The two solid lines above show the NLO and NNLO expansions, respectively, which show only a slightly weaker dependence on µ R . This feature was also seen in expansions of the NLL resummed cross section of Ref. [15], and should not be surprising, given that both the O(α s ) and O(α 2 s ) corrections have such a large size. Although the NNLO contribution includes logarithms of p T /µ R that help decrease the scale-dependence of NLO, the full NNLO contains additional terms whose scale dependence is compensated only at N 3 LO and beyond. The crosses again show the full NLO result, which is again in remarkable agreement with the NLO expansion of the resummed result. The upper In view of the large scale dependence still found at NNLO, one may wonder whether resummed perturbation theory may ultimately help to stabilize the predicted cross section with respect to scale changes. As described above (see discussion after Eq. (124)), the "minimal" Mellin inversion of the resummed cross section introduces unphysical contributions at z > 1 in Ω ab→c,resum (η, z) in (12), due to the presence of the Landau pole. In the present case, these contributions even turn out not to be controllable numerically. As remarked earlier, we leave for future work the implementation of a practical resummation formalism that includes full Mellin moment dependence, without unphysical support in z. Nonetheless, we can obtain finite and well-defined results by restricting z to the physical regime z < 1 in Ω ab→c,resum (η, z) and then in the convolution in Eq. (124). The upper lines in each of the two plots in the first row of Fig. 5 show the all-order result obtained in this way. As one can see, both of them are much flatter, showing very little residual renormalization scale dependence, and much reduced factorization scale dependence. This is precisely as would be expected from a full resummation formalism. We note that we find the same feature for RHIC energy at p T = 10 GeV. Our findings provide confidence that, once the unphysical regime z > 1 is adequately treated, resummation will yield valuable physical results. We caution that for the reasons just described the two results shown in the figures are not to be regarded as truly meaningful predictions of the resummation formalism we have developed here.
We finally examine the scale dependence for equal renormalization and factorization scales, µ R = µ F ≡ µ, as is often done in phenomenological studies. The lower left part of Fig. 5 shows  again the invariant cross section for E706 conditions, varying µ. The main patterns are as in the previous two figures. This time, we follow the results all the way down to µ = p T /2, although we do not necessarily favor such a small scale for the inclusive-hadron cross section: Given that the hadron takes only a fraction z ∼ 0.5 or so of the momentum of its parton progenitor, for a given hadron p T the hard-scattering typically will reside at a hard scale twice that value or so. The scale p T /2 may thus not reflect the hardness of the partonic interaction very realistically. Nonetheless, it is interesting to see that the various results edge closer together when the scale is chosen to have a small value. This becomes especially evident when going to RHIC energy in the lower right plot in Fig. 5. There, LO and NLO even meet at scale p T /2 for the value of p T = 12 GeV we consider, NNLO is only moderately higher, and the resummed result with z < 1 is also very close. It is interesting to note that such a tendency for the scale variation to narrow at low scales was observed in the literature also in various different contexts, for example early on in studies of prompt-photon production [74], but also recently for tt production at the LHC [75]. We stress again that we do not assign much significance to the precise location of the solid curve for the z < 1 resummed result. The fact that at low scales already the NNLO expansion is higher just shows once more that the result with z < 1 should be regarded as only a part of a fully resummed phenomenological cross section.

Conclusions
We have developed a threshold resummation for single-particle inclusive cross sections in hadronhadron scattering at next-to-next-to-leading logarithm, up to the ideal matching with exact nextto-next-to-leading order hard scattering functions. As in previous work on this subject, threshold resummation organizes leading-power plus distributions in the variableŝ 4 , the invariant mass of all radiation recoiling against the fragmenting parton.
New results include the definition and one-loop calculation of the matrix that organizes coherent soft radiation for all single-particle inclusive two-to-two partonic processes. This enables us for the first time to separate exact one-loop corrections to these processes between short-distance and long-distance factors, which are expanded in terms of the running coupling at hard and soft scales, respectively. The one-loop expansions of the factorized jet and soft functions that we derive reproduce all leading-power singular terms in the exact NLO calculations for partonic subprocesses. For completeness, the NLO singular terms are provided for each subprocess in appendices. In a test at phenomenologically-relevant kinematics, the O(α s ) expansions reproduce full NLO results to the accuracy of a few percent. The resummation analysis of this paper is given in Mellin moment space, and we compared with the closely-related NNLL resummed prompt-photon cross section of Ref. [18], developed using soft-collinear effective theory.
We have seen that the resummed single-particle inclusive cross section has a number of unique features that distinguish it from dihadron and single-photon inclusive cross sections. In particular, an inverse transform using the minimal prescription leads to unphysical contributions from z > 1 that we find are not numerically stable for single-particle inclusive cross sections, due to an enhanced unphysical range of the partonic fractional momenta for this process, which can cause partonic rapidities to become very large. We have provided an analytic Mellin inverse for Nindependent infrared scales, but given the apparent importance of threshold logarithms including Figure 5: Upper left: µ R -dependence of the invariant cross section for E706 energy at p T = 5 GeV for fixed µ F = 2p T . Upper right: µ F -dependence of the invariant cross section for E706 energy at p T = 5 GeV for fixed µ R = 2p T . Lower left: scale dependence of the invariant cross section for E706 energy at p T = 5 GeV, setting µ R = µ F ≡ µ. Lower right: scale dependence of the invariant cross section for RHIC energy at p T = 12 GeV, setting µ R = µ F ≡ µ.
N -dependence in the examples we studied, it seems natural to investigate the conventional "direct-QCD" approach further. Phenomenological applications thus will require further development, especially regarding the inverse transforms. This will be the subject of forthcoming work. We anticipate that the analysis given in this paper will be relevant to existing data at fixed-target energies, and to data from present and future hadronic colliders, in addition to the analysis of resolved photons at electron-hadron colliders. We expect that the formalism presented here will also have valuable applications in resummation studies of single-inclusive jet cross sections at hadron colliders.

Acknowledgements
We are grateful to Andreas Vogt for helpful communications and to Marco Stratmann for information on the DSS14 fragmentation functions, and to Andrea Ferroglia for a useful discussion. GS thanks the Pauli Center and the Institute for Theoretical Physics, ETH Zurich, as well as the Department of Physics, University of Vienna for their hospitality. The work of GS was supported in part by the National Science Foundation, grants PHY-1316617 and 1620628. FR is supported by the Department of Energy under Contract No. DE-AC0205CH11231, and the LDRD Program of Lawrence Berkeley National Laboratory. The numerical calculations presented here were in part carried out on the HPC resource bwUniCluster funded by the state of Baden-Württemberg.

A Anomalous dimensions
We present here the explicit low-order expansions of the various anomalous dimensions used for the resummed jet functions to the extent that they have not yet been given in the text. The functions A i ,B i andD i appearing in Eqs. (46), (70) and (71) are expanded as series in α s , The coefficients of A i are familiar. To NNLL, we use [56,76,77,78,79] A with N f is the number of flavors and The coefficientD (2) i was already given in Eq. (48). As shown in (50), it is directly related to the also widely used coefficient D (2) i . For completeness, we recall its explicit value [56,80]: In addition, for the recoiling jet in Eq. (71) we also need where as before b 0 = (11C A − 2N f )/(12π). The two-loop coefficientsB (2) i were already given in (63). As discussed in Eq. (64), they are obtained from the customary coefficients B (2) i which have been computed in [52]: Finally, we also present the expansion of the δ function contributions to the diagonal DGLAP splitting functions. These appear in Eq. (67) and also determine the factorization scale dependence of the hard function H in Eq. (37). Writing we have P and [81]

B Explicit NNLL forms of the exponents
Evaluating the exponent in Eq. (46) and choosing µ rf = √ŝ , one obtains an explicit expression for the NNLL expansion of the functionJ (i) in : where and Here b 0 , b 1 , b 2 are the first three coefficients of the QCD beta function which are given by [82,83] Note that we have obtained Eq. (142) by expanding the running coupling in the integrand of Eq. (46) as [56] α s (µ) = α s (µ R ) In this way, our perturbative expansion of the resummed exponents, which necessarily truncates the perturbative series, introduces dependence on a renormalization scale µ R ; see discussion after Eq. (115).
The corresponding result for the outgoing recoil jet may be written in compact form as where . (147)

C Ingredients for soft matrices
Every one-loop soft function for a given partonic channel is given by (see Eq. (112)) where in the last term the positive sign applies to all processes with a qq ( ) initial or final state, and the negative sign to all others. Below we present the matrices S (0) 0 and R 12 for all partonic channels. We also recall the one-loop soft anomalous dimension matrices Γ ab→cr,(1) (η), which may be found in this form in Ref. [27].
For qq → qq and qq → qq scattering we have The soft anomalous dimension matrix has already been given in Eq. (106). For qq → qq , qq → q q , and qq → qq scattering we have , Γ qq→qq, (1) where T = ln(1 − v) + iπ, U = ln(v) + iπ. For qq → gg we have For gg → qq we have the same S (0) 0 and Γ (1) as in (151), but For qg → qg and qg → gq we have the same S (0) 0 as in (151), but Finally, for gg → gg all three matrices have the block structure where, setting C A = 3 for simplicity,

D NLO expansions for other partonic channels
As before, we define L = ln(v) , For an arbitrary process ab → cr the first-order term in the product of the jet functions may be written as where C q = C F , C g = C A , and γ q = 3 2 C F , γ g = 2πb 0 = 1 6 (11C A − 2N f ) , The one-loop hard functions H (1) ab→cr used below may be found in Refs. [21,41,42].

D.1 qq → qq
The term Tr H (0) S (0) 0 is identical to that for the process qq → qq given in Eq. (119). Furthermore, and Tr with the trace terms for qq → qq also given in (119).

D.2 qq → qq
We have and In the last expression we have set C F = 4/3 and C A = 3 for simplicity.

D.3 qq → q q
We have Tr Tr H (0) Γ (1) † S In the last expression we have set C F = 4/3 and C A = 3 for simplicity.

D.4 qq → qq
We have Tr and Tr In the last expression we have set C F = 4/3 and C A = 3 for simplicity.

D.5 qq → gg
We have and Tr H (1) S (0) In the last expression we have set C F = 4/3 and C A = 3 for simplicity.
In the last expression we have set C F = 4/3 and C A = 3 for simplicity.