{'fraction_non_alphanumeric': 0.04133480825958702, 'fraction_numerical': 0.01801253687315634, 'mean_word_length': 4.2164839392190805, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sudakov Resummation and Finite Order Expansions of Heavy Quark Hadroproduction Cross Sections

We resum Sudakov threshold enhancements in heavy quark hadroproduction for single-heavy quark inclusive and pair-inclusive kinematics. We expand these resummed results and derive analytical finite-order cross sections through next-to-next-to-leading order. This involves the construction of next-to-leading order matching conditions in color space. For the scale dependent terms we derive exact results using renormalization group methods. We study the effects of scale variations, scheme and kinematics choice on the partonic and hadronic cross sections, and provide estimates for top and bottom quark production cross sections.


Introduction
Long-and short-distance dynamics in inclusive hadronic hard-scattering cross sections are factorized in Quantum Chromodynamics (QCD) into universal, non-perturbative parton distribution functions and fragmentation functions, and perturbatively calculable hard scattering functions. Remnants of long-distance dynamics in a hard scattering function can, however, become large in regions of phase space near partonic threshold and dominate higher order corrections. Such Sudakov corrections assume the form of distributions that are singular at partonic threshold. Threshold resummation organizes these double-logarithmic corrections to all orders, thereby extending the predictive power of QCD to these phase space regions.
Resummed cross sections constitute an approximate sum of the complete perturbative expansion if, at each order, the Sudakov corrections dominate. The numerical evaluation of resummed cross sections requires a prescription to handle infrared renormalon singularities. Resummed cross sections may also be expanded to provide estimates of finite higher order corrections which do not suffer from renormalon problems. In this paper we shall employ the resummed cross sections in the latter fashion: as generating functionals of approximate perturbation theory.
Resummed results for heavy quark production at leading logarithmic accuracy have been presented some time ago [24,25,26,27,28,29]. In our paper we expand the NLL resummed cross sections presented in Refs. [4,5,12] and derive complete analytic expressions through next-to-next-to-leading order (NNLO) for double-differential heavy quark cross sections in two different kinematics: heavy quark pair-inclusive and single-heavy quark inclusive. To achieve next-to-next-to-leading logarithmic (NNLL) accuracy we include color-coherence effects and contributions due to soft radiation from one-loop virtual graphs via matching conditions. Our paper is organized as follows. In section 2 we discuss both types of kinematics and their singular functions. Section 3 describes the construction of the resummed cross sections. In section 4 we expand the resummed cross section to NLO and NNLO and present analytical NNLL double-differential cross sections at each order, with both types of kinematics, in the gluon-gluon (gg) and quark-antiquark (qq) channels. We numerically study the inclusive partonic and hadronic cross sections in sections 5 and 6 respectively. Our conclusions are presented in section 7. Appendix A contains the NLO matching terms for pair-inclusive kinematics. In appendix B we collect all our explicit expressions for the NLO and NNLO differential cross sections.
A companion study [30], also addressing third and fourth order corrections, contributions from subleading logarithms, and some differential distributions, already featuring some of the second order results for inclusive cross sections derived in this paper, was recently presented by one of us.

Kinematics and cross sections
The kinematics of inclusive heavy quark hadroproduction depend on which final state momenta are reconstructed. In threshold resummation this kinematics choice manifests itself at next-toleading logarithmic level [12]. We discuss two types of near-elastic kinematics in heavy quark hadroproduction, one-particle inclusive (1PI) and pair-invariant mass (PIM) kinematics.

One-particle inclusive (1PI) kinematics
Heavy quark hadroproduction in 1PI kinematics is defined by where h 1 and h 2 are hadrons, X[Q] denotes any allowed hadronic final state containing at least the heavy antiquark, and Q(p 1 ) is the identified heavy quark with mass m. The hadronic invariants in this reaction are S = (P 1 + P 2 ) 2 , T 1 = (P 2 − p 1 ) 2 − m 2 , U 1 = (P 1 − p 1 ) 2 − m 2 , and where S 4 is a measure of the inelasticity of the hadronic reaction (1). Near threshold, reaction (1) is dominated by the partonic subprocesses If X ′ [Q](p ′ 2 ) = Q(p 2 ), the reaction is at partonic threshold and the heavy antiquark has momentump 2 . Note that threshold production does not mean that the heavy quarks are produced at rest. The qg and qg channels contribute at one order higher in α s than the reactions (4) and (5) 2 . The partonic invariants corresponding to (2) are The invariant s 4 = (p ′ 2 ) 2 − m 2 which measures the inelasticity of the partonic reactions (4) and (5) is related to the other partonic invariants by s 4 = s + t 1 + u 1 .
The inclusive partonic cross section may be calculated from where β = 1 − 4m 2 /s and The recoil momentum p ′ 2 may be split into the momentum at threshold,p 2 , and the momentum of any additional radiation above threshold, k, i.e. p ′ 2 =p 2 + k. Then, when k 2 is small, we can define a dimensionless weight w 1PI [3] that measures the distance from threshold in 1PI kinematics which can in turn be expressed in terms of a vector ζ µ 1PI : Ref. [33] contains an exact NLO treatment of this kinematics at the parton and hadron levels.
At the parton level, the important reactions are q(k 1 ) + q(k 2 ) −→ QQ(p ′ ) + X ′ (k) , where p ′2 = M 2 is the pair-mass squared. If X ′ (k) = 0, the reaction is at partonic threshold with M 2 = s. Then where β M = 1 − 4m 2 /M 2 and θ is the scattering angle in the parton center-of-mass frame.
In PIM kinematics, the inclusive partonic cross section may be calculated from with The weight w PIM measures the inelasticity in PIM kinematics, where z ≡ M 2 /s. Equation (17) defines the vector ζ µ PIM in terms of the heavy quark pair momentum at threshold (indicated by the bar). An exact NLO treatment of this kinematics at the parton and hadron levels may be found in Ref. [34].

Singular functions
The double-differential hadronic cross section for reactions (1) and (11), d 2 σ h 1 h 2 /dV dW , enjoys the factorization [35] where ω ij is the partonic cross section, or hard scattering function, whose dependence on kinematics is indicated by the weight w K with K = (1PI, PIM). The variables V and W represent either 1PI or PIM kinematic variables such as the transverse momentum and rapidity of either a single heavy quark or the heavy quark pair, respectively. The parton momentum fractions x 1 and x 2 have lower limits x − 1 and x − 2 which depend 3 on V and W . We shall resum the higher order logarithmic contributions to ω ij that are singular at threshold. The arguments of these logarithms are the weights w K . Thus, the 1PI singular functions are plus distributions in s 4 , while the PIM singular functions are plus distributions in 1 − z These singular functions yield finite large results when convoluted with smooth but rapidly changing functions such as parton densities. Note that we have normalized the 1PI functions to have mass dimension −2. Because our paper deals mostly with finite order cross sections, we denote corrections as leading logarithmic (LL) if l = 2i + 1 at order O(α i+3 s ), i = 0, 1, . . ., as next-to-leading logarithmic (NLL) if l = 2i, etc.
It is often convenient to work in moment space, defined by the Laplace transform with respect to w Kf The upper limit of this integral is not very important, and may be set to 1, where ln w K = 0. Under Laplace transformations, the plus distributions in Eqs. (22) and (23) become linear combinations of ln k (Ñ) with k ≤ l + 1 andÑ = N exp(γ E ) where γ E is the Euler constant 4 . The precise correspondence to second order can be found in Ref. [2], and through fourth order in Ref. [30]. We shall almost always talk about logarithmic accuracy in momentum space. Occasionally in the following we denote in moment space leading logarithmic corrections at order O(α i+2 s ), i = 1, 2, . . . to correspond to l = 2i, NLL ones to l = 2i − 1 etc. Although the Laplace transformation of the l'th plus distribution generates lower powers of ln N besides ln l+1 N, there should be no confusion in practice.
We work in axial gauge, n · A = 0, with the gauge vector chosen as n µ = ζ µ K . The implicit renormalization scheme is that of Ref. [36] and m is a pole mass. The renormalization scale µ R is assumed to be equal to the factorization scale µ, except where explicitly indicated otherwise. The scale dependence of the coupling constant is controlled by the QCD β-function S where Y is the pair rapidity. 4 All the N 's in the remainder of this work are actuallyÑ , unless specified otherwise.

Refactorization
The resummation of the singular functions in Eqs. (22) and (23) in the perturbative expansion of ω ij rests upon the refactorization of ω ij into separate functions of the jet-like, soft, and off-shell quanta that contribute to its quantum corrections. This refactorization, valid in the threshold region of phase space, is pictured in Fig. 1. Each of the functions ψ, h, and S organizes Figure 1: Refactorized form of heavy quark partonic cross section near threshold.
large corrections corresponding to a particular region of phase space. The meaning of the I, J indices will be given shortly and the coefficients c ij I are given below in Eq. (43). Factorizations of this type have been discussed earlier for deep-inelastic scattering and Drell-Yan production [1], heavy quark hadro-and electroproduction [4,5,8,13], dijet [6] and Higgs [16] production, and single-particle/jet inclusive cross sections [12,21]. They may be generalized to include recoil effects [37]. Figure 1 indicates that each factorized function depends on its own weight function [3]. In essence, the choice of working in 1PI or PIM kinematics depends on how the total weight w K is constructed from the individual contributions of each weight function [6,12]: 1PI : In the case of 1PI kinematics it is convenient to define 1PI : In both kinematics, the moments of the ij partonic cross section, Eqs. (4), (5), (12), and (13), can be written in refactorized form, up to O(1/N) corrections, as [4,5,6,12] ω ij (N, s, The indices 5 I and J take values in the color-tensor space spanned by the invariant tensors that combine the SU(3) representations of the external partons at threshold into a singlet. The vector ζ µ K defines the kinematics, see section 2. The "incoming-jet functions"ψ i/i describe the dynamics of partons moving collinearly to the incoming parton i. The distributionsφ i/i are defined at fixed light-cone momentum while the functionsψ i/i are defined at fixed k i · ζ K . The (real)ψ i/i functions include all leading and some next-to-leading singular functions and are diagonal in color-tensor space. The functionS ij JI , a Hermitian matrix in color-tensor space, summarizes the dynamics of soft gluons that are not collinear to the incoming partons and contains the remaining next-to-leading contributions. Note that (next-to-leading) singularities associated with soft radiation from the outgoing heavy quarks 6 are included inS ij JI [4,5]. Finally, the Hermitian matrix H ijIJ ≡ h ijI h * ijJ incorporates the effects of far off-shell partons and contains no singular functions.
The jet and soft functions in Eq. (29) can each be represented as operator matrix elements [4,5,6,7,12]. The refactorization of the cross section in Eq. (29) can in fact be seen as a separation of near-threshold degrees of freedom into distinct effective field theories [3,4,5]. In each effective theory a resummation of the singular functions may be performed via appropriate evolution equations. The resummed cross section is composed of contributions from all these effective theories which must be matched together properly, at a specified scale and to a certain order in perturbation theory. In this paper we match to NLO, leading to NNLL accuracy in our finite order expansions. Our matching scale is the heavy quark mass m. The matching procedure is described in section 4.
In the following we describe the resummation of the jet and soft functions. Most of these steps have been discussed previously in the literature [4,5,6,12] so that our description is brief.

Resummed Jet and Soft Functions
The (exponentiated) N dependence of the ratioψ i/i /φ i/i in Eq. (29) follows from its factorization properties [1,3]. The functionψ i/i obeys, beside a renormalization group equation, an evolution equation governing the energy dependence, expressed as gauge-dependence [1,3,38,39]. Solving these equations leads to the resummed form of this ratio. The function ψ has been defined in PIM kinematics for i = q [1] and g [6] with the 1PI equivalents given in Refs. [13,12]. The resummation of Sudakov logarithms in 1PI kinematics was verified [13] to trace those in PIM kinematics [1]. Expressed in terms of w i , theψ i/i are the same in PIM and 1PI kinematics.
with the MS exponent The functions A i and ν i differ by color factors for i = q and g so that with K = C A (67/18 − π 2 /6) − 5n f /9 [40].
To incorporate the effects of scale changes in finite-order α s expansions, we note that the ratioψ i/i /φ i/i transforms under renormalization scale, µ R , and factorization scale, µ, as where γ i and 2γ i/i are, respectively, the anomalous dimension of the quantum field i and of the operator whose matrix element represents [41,42] the MS-densityφ i/i . The anomalous dimension for the renormalization scale dependence ofψ i/i [6] in Eq. (33) does not depend on the moment N because the ultraviolet (UV) divergences ofψ i/i are only due to wavefunction renormalization of parton i. Hence γ ψ i = 2γ i where γ i is calculated in the axial gauge. The anomalous dimension γ i/i controlling the factorization scale dependence does depend on N.
The axial gauge anomalous dimensions are (neglecting O(1/N) terms) where b 2 is given below Eq. (25) and K is given below Eq. (32).
The composite operator [7] that definesS ij JI contains UV divergences beyond the self energies of its external legs. An additional renormalization removes these extra divergences. However, the product H ijIJSijJI has no UV divergencies beyond those taken into account by the renormalization of the elementary fields and couplings of the theory, because the renormalization-group invariant cross section and this product only differ by the functionsψ i/i which are renormalized as if they were (the "square" of) elementary fields (33). Hence the extra UV divergencies ofS ijJI are balanced by similar ones in H ijIJ . Therefore, the unrenormalized expressionsS bare ijJI and H bare ijIJ renormalize multiplicatively [4,5,6] The factors Z i denote the renormalization constants of the external fields i = q,q, g. For a given parton pair ij, the (Z ij S ) IJ constitute a matrix, in color-tensor space, of renormalization constants for the overall renormalization of the soft function. Note that (Z ij S ) IJ also includes the wave function renormalization of the external heavy quark legs. Equation (40) leads directly to a renormalization group equation for the matrixS ijJI [4,5,6] The soft anomalous dimension matrix Γ ij S,IJ is obtained from the renormalization constants (Z S ) IJ , computed in d = 4 − ǫ dimensions, by [43] with g 2 s = 4πα s . The solution of Eq. (41) is in general expressed in terms of path-ordered exponentials [4,5,6], see below.
An explicit expression for Γ ij S,IJ requires a choice of basis tensors in color-tensor space. We choose an s-channel singlet-octet basis where T c F are the SU(3) generators in the fundamental representation, where the indices m, n, k, l take values. The indices a, b, c are adjoint indices while d abc and f abc are the totally symmetric and antisymmetric SU(3) invariant tensors, respectively. In this basis, the one-loop soft anomalous dimension matrix for the qq process is [4,5,22] with matrix elements The function L β is defined as The variables ν i , i = 1, 2 are The one-loop soft anomalous dimension matrix for the gg process [5,22] can be written as where the three independent matrix elements are Note that we did not absorb the function ν i of Eq. (31) into Γ S . Combining the solution of Eq. (41) with the resummed jet functions and the hard function, and using matrix notation for H ij andS ij , the full resummed partonic cross section is where the trace is in color-tensor space, and P refers to path-ordering in µ ′ such that Γ ij S (α s (m)) is ordered to the far right while Γ ij S (α s (m/N)) is ordered to the far left. The operationP orders in the opposite way. The exponential exp(Ẽ i ) is given by The exponential exp(Ẽ j ) is identical except for obvious relabelling. In Eq. (50) we used the renormalization group behaviour, derived from Eqs. (39) and (40), of the product H ij IJSij JI where we have only indicated the µ R and N dependence. Notice that all µ R dependence cancels in Eq. (50). So far, the results have all been given in the MS-scheme. The qq results can also be presented in the DIS scheme, in which the functionẼ q of Eq. (31) is In some cases, e.g. for the unexpanded resummed cross section, it is preferable to diagonalize the matrix Γ ij S,IJ [7,11,22,44] by a change of basis in color-tensor space, written in matrix notation with Γ R,S diagonal. The channel labels have been suppressed. The columns of the matrix R are the eigenvectors of Γ S . This leads to ordinary exponentials in the solution of Eq. (41), rather than path-ordered ones. The elements on the diagonal are the (complex) eigenvalues λ I . In the new basis, indicated by the subscript R, the matricesS and H areS again in matrix notation. Eq. (41) then becomes The matrix R qq is given in Refs. [11,22] while R gg appears in Ref. [22]. The solution of Eq. (57) isS so that (reinstating the channel labels) ω res ij (N, s, t 1 , u 1 , m 2 , µ 2 , α s (µ)) = The exponential in Eq. (59) carrying color-tensor indices is with exp(Ẽ i ) and exp(Ẽ j ) given in Eq. (51). We derived the bulk of our subsequent results in both bases.

NNLL-NNLO expansions for partonic cross sections
In this section we derive analytical expressions for partonic double-differential heavy quark cross sections up to NNLO by expanding the resummed cross section of the previous section. We concentrate here on the derivation; our explicit results are collected in the appendices. We obtain expressions for both 1PI and PIM kinematics and in both MS and DIS factorization schemes. Our aim is NNLL accuracy, as defined in section 2.4, including the scale dependent sector where we distinguish between the renormalization scale µ R and mass factorization scale µ. This means that for coefficients of ln i (µ 2 /m 2 ) or ln i (µ 2 R /m 2 ) we include the most singular plus distribution and the two next-most-singular ones.
To achieve NNLL accuracy we must derive NLO matching conditions for certain functions in the resummed cross section, Eq. (50). To this end we expand the factors in this equation and identify order by order the perturbative coefficients of the functions in Eq. (50). The coefficients can in part be explicitly computed and in part inferred by matching to exact NLO results.

Expansions
Let us expand each factor in Eq. (50) in powers of α s . We set µ R = m here because Eq. (50) is independent of µ R , leading to a cross section expansion in α s (m) which can easily be changed back to α s (µ R ). Neglecting 1/N terms, the two-loop expansion of exp(Ẽ i ) may be written as The coefficients C i,(n) k can be computed from the results in section 3. To NNLL accuracy they are where we have used the notation and where A (n) i and ν i,(n) are the coefficients of (α s /π) n in the expansion of the functions A i and ν i in Eq. (32). The two-loop expansion of the path-ordered exponential in Eq. (50) reads P exp is the coefficient of α s /π in the one-loop soft anomalous dimensions of Eqs. (44) and (48). TheP ordered exponential is identical to Let us for the moment choose PIM kinematics by putting N t = N u = N. We substitute Eqs. (61), (64), and (65) into Eq. (50), finding the moment space expression for the NNLL-NNLO expanded cross section: With 1PI kinematics we use Eq. (28) and expand ln(N u ) = ln(N) + ln(−u 1 /m 2 ) etc. To transform from moment to momentum space we use results given in the appendices of Refs. [2,13,30]. Note that we do not keep subleading terms in our expansions when we invert to momentum space [30].

Matching
We can identify the functions we must determine from Eq. (66). The coefficients C i,(n) k are given in Eq. (62). The matrices Γ ij (1) IJ are given in Eqs. (45) and (49). The remaining unknowns are The leading-order soft functionS where {c I } is our basis in color-tensor space. Recall that we have chosen the s-channel singletoctet basis, Eq. (43). Note thatS (0) = S (0) . The lowest-order hard function H is calculated by projecting the lowest order amplitude A ij onto this color basis: where the trace acts in ordinary color space. Note that at this stage the function h ij I may still depend on all other quantum numbers of the process under consideration, such as Lorentzspinor and spin degrees of freedom. We are however not interested here in spin-dependent observables and we trace over all remaining degrees of freedom. The matrix H (0) ij IJ is real and symmetric but may have zero eigenvalues if the Born process does not span the full color-tensor space, as is the case for the qq channel. For the qq channel [11] the results arẽ with H (0) qq 22 given by The lowest-order results are kinematics-independent. For the gg channel we find where and It now remains to determine 7 Tr H We do this by matching. The necessary matching conditions can be inferred by comparing the expansion in Eq. (66) to exact results for the partonic cross section. At lowest order, the matching condition is It is straightforward to check that this condition is fulfilled by our results forS ij IJ . Neglecting 1/N terms, and choosing the case K = PIM for the moment, we have at NLO the matching condition Let us discuss this equation. To begin, the first term on the right hand side may be determined from the exact results in Refs. [32,33,34], converted to moment space. To NNLL accuracy this term consists of all terms in the differential one-loop partonic cross sections that either diverge logarithmically or are constant. The constant terms consist, first of all, of the virtual graph contributions. In addition there is the soft-gluon contribution from the radiative graphs. The virtual contribution is kinematics independent but the soft gluon one is not because what is defined as soft differs for 1PI and PIM kinematics. The respective criteria are s 4 < ∆ and 1 − z < δ, with ∆ and δ defined in Eqs. (22) and (23). The virtual graphs contain various divergences. The ultraviolet ones are cancelled by renormalization of the QCD coupling constant and the heavy quark mass, for both kinematics in the scheme of Ref. [36]. Infrared divergences are cancelled by adding the virtual and soft contributions, so that only collinear divergences remain. These must be subtracted via mass factorization. The corresponding subtraction terms in general involve convolutions of Altarelli-Parisi splitting functions [45] with the lowest order cross section. However, we must mass factorize only with the soft plus virtual parts of these splitting functions. These depend on kinematics, see e.g. Eqs. (6.8) and (6.13) in [32]. What remains at NLO after mass factorization can now be categorized as terms multiplying ln i (∆/m 2 ) or ln i (δ), i = 1, 2, and other terms. The logarithms ln(∆/m 2 ) or ln(δ) should be completed to plus-distributions via Eqs. (22) and (23). The final result of these procedures, in moment space, is what constitutes the first term on the right hand side of the equals sign in Eq. (80). The other terms on the right hand side merely subtract all terms containing the singular functions in Eq. (22) or (23).
For 1PI kinematics there is a slight subtlety, related to the use of N t and N u in Eq. (28). Expanding ln(N t ) = ln(N) + ln(−t 1 /m 2 ) (and ln(N u ) likewise) leads to extra terms containing ln(−t 1 /m 2 ) and ln(−u 1 /m 2 ) but not ln(N) and one must be careful in accounting for such terms.
We can now obtain specific results for the matching terms in momentum space. At NLO, T (1) ij K on the left hand side of Eq. (80) for 1PI kinematics, is the expression for s 2 d 2 σ (1),S+V ij (s, t 1 , u 1 )/dt 1 du 1 given by Eq. (4.7) of Ref. [33] for ij = qq, and Eq. (6.19) of Ref. [32] for ij = gg, modified by removing a factor α s /π as well as all terms containing ln 2 (∆/m 2 ) or ln(∆/m 2 ) and terms containing ln(µ/m) or ln(µ R /m) 8 . The definition of the hatted symbol in Eq. (81) makes the delta function explicit and is useful when we present the two-loop results.
In PIM kinematics the NLO matching terms may be determined as follows. We start from the above 1PI results and modify those parts that are kinematics dependent. From our discussion following Eq. (80), it is clear that these are only the soft contributions and the mass factorization subtraction terms. We shall label these terms by the superscript S+MF. The soft contributions to the PIM NLO cross sections can be computed using the results of Ref. [34]. As mentioned in the above, we mass-factorize these results with only the soft plus virtual part of the splitting functions in PIM kinematics. The same procedure is done for the 1PI case using the results of Refs. [32,33]. For channel ij we then compute the MS quantity where the subscript PIM on the second term on the right hand side indicates that for t 1 , u 1 one must use the expressions in Eq. (14). Note that we may replace β M by β to the accuracy at which we are working in this paper, and we shall do so for all following results. Explicit expressions for the terms on the right hand side of Eq. (82) can be found in appendix A. Then the PIM equivalent of Eq. (81) is, accounting for an overall jacobian factor originating from the definition of the cross section, Note that terms containing renormalization and factorization scale dependence can be easily included by expanding the corresponding scale dependent exponentials in Eq. (50) and adjusting the coefficients in Eq. (66) correspondingly. Alternatively, the scale dependence can be constructed using renormalization group methods, which we do in the next section.
Changing to the DIS scheme at NNLL for the qq channel involves using Eq. (53) rather than Eq. (31) in the expansion, Eq. (61), and extra constant terms at NLO as given in appendix B.
Because of their length, we have collected all our results in appendix B. In the next section we perform a numerical study of the results obtained in this section for the inclusive partonic cross sections.

Scaling functions
It is convenient to express the inclusive partonic cross sections in terms of dimensionless scaling functions f (k,l) ij that depend only on η, defined in Eq. (20), as We derive LL, NLL, and NNLL approximations to f (k,l) ij (η) for both the qq and gg channels using the results of the previous section, via Eqs. (8) and (15). Where possible we will compare with exact results. We also examine the results for the qq channel in the DIS scheme. Of course, a complete comparison between the DIS and MS schemes also requires the use of DIS parton densities, as discussed in the next section.
We begin by showing the scaling functions f (k,0) ij (η), k = 0, 1, 2, in both 1PI and PIM kinematics. In Fig. 2 we present the qq Born and MS NLO results for the scaling functions f  qq , the LL approximations already reproduce the exact curve quite well in both kinematics [24], but the NLL approximations agree much better with the exact results to larger η. Adding the NNLL corrections, numerically dominated by the negative contributions from gluon exchange between final state heavy quarks (Coulomb terms), leads to excellent agreement with the exact curve over a large range of η. Such progressive improvement was also observed, for 1PI kinematics, in Ref. [46] where the absolute threshold limit was taken. We also see that, while the LL approximations in 1PI and PIM kinematics differ and under-or overestimate the true result, the NLL approximations in both kinematics agree better with each other, and even more so at NNLL, at least for η ∼ < 0.1. To summarize the one loop results in this channel, differences due to kinematics choice decrease as the logarithmic accuracy increases.
The same conclusions hold for the NLO scaling function in the gg channel, f (1,0) gg , shown in Fig. 3. The Coulomb terms are positive here. At large η, the agreement between the 1PI and PIM results is not as good as in the qq channel: the 1PI results at NNLL remain positive while the PIM ones are negative. gg for both kinematics. Again, the LL approximations differ but there is good agreement between 1PI and PIM kinematics at NLL accuracy. We also display the NNLL approximations for f (2,0) qq and f (2,0) gg , presently the best estimates for these functions, at least for η ∼ < 1, since exact calculations are not yet available. Near threshold, the NNLL approximations show some deviation from the NLL approximation due to interplay between the one-loop LL contributions and the one-loop Coulomb terms. Note that while the qq channel shows relatively good agreement between the two kinematics at large η, the gg results are quite different since the 1PI results are positive while the PIM results are negative.
We have also calculated the behavior of the scaling functions f (k,0) ij (η), k = 1, 2, in the threshold limit s → 4m 2 , keeping only terms that grow as ln(β) to next-to-leading accuracy. The resulting expressions are given in appendix B.5, and shown in Figs. 2(b), 3(b), 4(b) and 5(b). They agree fairly well with the NLL PIM results.
In general, comparison of the exact and approximate results in the various channels and kinematics in Figs. 2 and 3 indicates the range of η values over which we expect the threshold approximation to be valid. The large-η behaviour of the scaling functions, furthest away from threshold, may be approximated by the methods of Ref. [47,48,49,50]. However, the behavior We now check if the qq MS results are mirrored in the DIS scheme. Since the scaling functions in the two schemes differ already at LL, the higher-order scaling functions will be quite different numerically. Of course, this difference should be compensated in principle by the parton densities for the scheme-independent hadronic cross section. Figures 6 and 7 show the functions f We next discuss the functions controlling the scale dependence. They can be determined exactly at NLO and NNLO using renormalization group methods. The exact f . We neglect flavor mixing terms which are of order 1/N, where N is the moment variable. We checked that at Tevatron energies the error due to this approximation is less than 1% at NNLO. In this way we obtain where P (0) qq , P (0) gg and P (1) qq , P (1) gg are the one-and two-loop splitting functions [51,41,52], and ≃ indicates the neglect of flavor-mixing terms. The convolutions involving a scaling function f with η(x) = xs/(4m 2 ) − 1. The standard convolution of two splitting functions, P (0,0) qq and P (0,0) gg , in Eqs. (87) and (90) are with i = q, g. Equations (85) and (88)  gg in Ref. [31]. We have also checked that the above results agree to NNLL with the expressions in appendix B when integrated as in Eqs. (8) and (15). Note that in appendix B we have also given results for terms in the differential NNLO functions controlling the scale dependence beyond NNLL accuracy, thus deriving all the soft plus virtual terms in these functions. The exact hard terms are calculated only for the integrated cross section as above.
We begin by showing the scale-changing scaling functions in the qq channel and MS scheme, comparing 1PI and PIM kinematics. In Fig. 8 we show f (1,1) qq and note that the PIM LL approximation reproduces the exact curve somewhat better than the 1PI LL aproximation. The NLL approximations agree better, even for larger η.
In Fig. 9 we show the NNLO scaling functions f (2,1) qq and f (2,2) qq . We compare the exact curves calculated from Eqs. (86) and (87) with our LL, NLL, and NNLL approximations. Again we see that the NNLL approximations provide a remarkably good description of the exact results, both in shape and magnitude. The NNLL curves for 1PI and PIM kinematics are in very good agreement with each other, i.e. ambiguities from the kinematics choice are very mild. Similar conclusions hold for the gg scaling functions f (1,1) gg , f (2,1) gg , and f (2,2) gg , shown in Figs. 10 and 11. Note however that also here the NNLL 1PI and PIM results for f (2,1) gg differ at large η. The PIM NNLL scaling function differs significantly from the exact result.
For completeness, we also display the NNLO scaling functions f

Hadronic cross sections
In the previous section we examined the quality of our threshold approximations at the parton level. Here we assess these approximations at the hadron level for inclusive top quark production at the Fermilab Tevatron and bottom quark production at HERA-B. The inclusive hadronic cross section is the convolution, Eq. (19), of parton distribution functions with the partonic cross section, expressed in terms of the scaling functions, Eq. (84). To facilitate the understanding of the results in this section in terms of those of the previous section, we plot the flux factors Φ ij (η, µ 2 ), Eq. (18), for the above cases in Fig. 13. They show which η values receive the most weight in the convolution integral. In our numerical studies we use the two-loop expression of α s and the CTEQ5M (MS scheme) or CTEQ5D (DIS scheme) parametrizations of the parton distributions [53] not only for the NLO results, but also for the NNLO (NNLO parton distributions are not yet available) and LO results. Thus in this section we keep the nonperturbative part of our results fixed when studying the effect of increasing the perturbative order of our partonic cross sections. For top quark production at the Tevatron and bottom quark production at HERA-B the calculations probe the moderate to large x region where the parton distributions are well known, see Fig. 13. We use five and four active flavors respectively for these cases and fix µ R = µ. Except where specified otherwise we have multiplied the non-exact scaling functions at NLO and NNLO with a damping factor 1/ √ 1 + η, as in Ref. [54] in order to lessen the influence of the large η region of the scaling functions where threshold logarithms become less dominant and we lose some theoretical control 9 . Figure 14 demonstrates that this factor indeed damps the large η region of the NNLO scaling functions f (2,0) ij , while leaving the small and medium η regions unaffected. The effect of the damping factor will also be made explicit in the tables.

Results for tt production at the Tevatron
We first discuss top quark production at the Tevatron in proton-antiproton collisions. (For a review of the data see Refs. [55,56].) We give results for both Run I ( √ S = 1.8 TeV) and Run II ( √ S = 2.0 TeV). At the Tevatron, top quarks are mainly produced in pairs. With a top quark mass of 175 GeV, the dominant production channel is qq annihilation, constituting about 90% of the total cross section at the Born level for √ S = 1.8 TeV, with the gg channel making up the remainder. The NLO corrections in the qq channel are moderate, of the order of 20%, whereas those in the gg channel are more than 80% so that the gg channel gains significance at NLO 10 . These large corrections originate predominantly from the threshold region. Figure 13 shows that the range 0.1 ∼ < η ∼ < 2 contributes the most. At still higher orders, this trend continues: the relative corrections to the gg channel are larger than those for the qq channel, as we shall see. This is due to the larger color factors in the analytical expressions for the corrections for the gg channel. As we mentioned earlier, the qg andqg channels give negligible contributions to the total cross section and are not considered here.
We begin by comparing MS results in 1PI and PIM kinematics, including only the qq channel in Eq. (19), at √ S = 1.8 TeV. In Fig. 15 we show the Born cross section and the exact and approximate NLO corrections, the latter at both NLL and NNLL accuracy, for µ = m in the range 150 < m < 200 GeV. We see that our NLO 1PI approximations are a little larger than the exact result while the NLO-NNLL approximation in PIM kinematics is indistinguishable from the exact answer. The corrections are about 20-30% in the mass range shown. In Fig. 16 we give the equivalent results for the gg channel. Here the 1PI approximations agree with the exact result better than the PIM ones. Figure 17(a) displays the approximate qq NNLO corrections at NLL and NNLL accuracy for µ = m as a function of the top quark mass in a direct comparison of the 1PI and PIM results. The shaded area indicates the kinematics ambiguity at NNLL, about 0.6 pb at m = 175 GeV. The figure shows that the NNLL ambiguity is larger than the NLL one. (Thus the small size of the NLL kinematics ambiguity seems somewhat accidental.) Figure 17(b) shows similar results for the gg channel.
For completeness, we show the corresponding results for the qq channel in the DIS scheme in Figs. 18 and 19. Again PIM kinematics approximates the exact results somewhat better than 1PI kinematics at NLO-NNLL. We see that the NNLL-NNLO kinematics ambiguity in Fig. 19, again indicated by the shaded region, is greatly reduced compared to the MS case, Fig. 17(a). It should however be kept in mind that, in the DIS scheme, the parton densities absorb large threshold logarithms which are not properly accounted for at NNLO if one uses NLO parton distributions as in Fig. 19. Therefore it seems likely that the DIS scheme kinematics ambiguity is somewhat underestimated in Fig. 19.
In Fig. 20 the sum of the qq and gg channels in the MS scheme is shown as a function of the top quark mass. We display the exact NLO cross section and the approximate NNLO cross section, which is the sum of the exact NLO cross section and the NNLL-NNLO corrections. We show results for µ = m/2, m, and 2m. We see that the NNLO cross section is uniformly larger than the exact NLO one, although less so in PIM kinematics, and that the scale dependence of the NNLO cross section is considerably reduced relative to NLO. Comparing Figs. 20(a) and (b) we observe that the kinematics ambiguity is larger than the scale uncertainty.
We now turn to a more detailed study of the scale dependence of the inclusive MS top      cross section. We show the sum of the qq and gg channels in Fig. 21 at √ S = 1.8 TeV and m =175 GeV at several orders and to different accuracies. The Born and NLO results shown are exact. One of the NNLO curves is constructed by adding the contributions from the NNLL approximate two-loop scaling functions f (2,l) ij , l = 0, 1, 2, to the exact NLO results, the other by adding instead the exact f (2,1) ij and f (2,2) ij functions to the approximate f (2,0) ij , thus making the changes to the partonic cross section when changing µ exact. The differences between these two NNLO curves are due to subleading terms and represent, for each kinematics, the corresponding ambiguity in the scale dependence. At very small µ the contributions of the terms involving f (2,1) ij and f (2,2) ij are much larger that the contribution from f (2,0) ij . The sizable difference between the two NNLO curves in Fig. 21(a) is in fact mainly due to the difference between the exact and NNLL 1PI results for f (2,1) qq at medium and large η in Fig. 9(a). Even if we include all soft plus virtual terms in the approximate 1PI f (2,1) qq , as derived in appendix B, there is still a sizeable difference from the exact result. Therefore this difference stems mainly from the hard, i.e. O(1/N), terms in f (2,1) qq . The NNLO differences in PIM kinematics are much smaller than in 1PI kinematics, in correspondence with the good agreement of the exact and NNLL PIM results for f (2,1) qq over all relevant η, as shown in Fig. 21(b). Next, in view of the upgrade in energy for the Tevatron from √ S = 1.8 to 2.0 TeV, we investigate the √ S dependence of the top quark production cross section. In Fig. 22 we present the inclusive cross section for the sum of the qq and gg channels in the MS scheme at NLO and NNLO as a function of √ S. The NNLO curves result from adding the NNLL NNLO corrections to the exact NLO cross section. We have normalized all calculations to the value of the exact NLO cross section at µ = m. Comparing 1PI and PIM kinematics we find that at lower energies, where the qq channel is dominant, the NNLO results are 10-30% larger than at NLO in both kinematics. As √ S increases, the gg channel, with its larger corrections, grows in importance. The cross sections are also more sensitive to the large η behavior of the scaling functions. This leads to large kinematics differences at energies above 5 TeV although the scale uncertainties remain small. In PIM kinematics this is due to f (2,0) gg which is large and negative for η > 0.1 at NNLL, thus reducing the NNLO results relative to the exact NLO cross section.  Fig. 21 where we also show NNLO results corresponding to the EX+DF choice (dash-dotted lines). Comparing the NNLO and NLO cross sections in the tables shows that the NNLO corrections are much smaller in PIM kinematics than in 1PI kinematics, as already shown in Figs. 22(a) and (b). We see that for the top cross section the effect of the damping factor is rather small. We also note that the use of exact f (2,1) ij and f (2,2) ij affects the µ = m/2 case the most. The difference between the exact and approximate calculation is larger for 1PI than PIM kinematics, as can also be observed in Fig. 21.
We note that AP results in 1PI kinematics were also presented in Ref. [30]. There are some small numerical differences with the results in this paper stemming mainly from using slightly different analytical expressions, all equivalent at threshold. To be specific, the expressions t 2 1 + u 2 1 and s 2 − 2t 1 u 1 are equivalent at threshold and either choice can be made in our NNLO expansions. Different choices produce slightly different numerical results away from threshold and the variation in these results represents a small but inherent uncertainty in the cross sections.
Let us finally comment on the applicability of our results for top quark production at the LHC where √ S = 14 TeV. At this pp collider the gg channel is dominant (about 90% of the   total) because only sea quarks contribute to the antiquark distributions in the qq channel.
Although top production at these energies is far from the hadronic threshold region, it might be close enough to partonic threshold for threshold resummation to be relevant since the gluon flux may favor small values of s = x 1 x 2 S ≪ S (see the discussion in Refs. [57,58]). Fig. 23 indeed confirms this, but also shows that the hadronic cross section is sensitive to the highenergy behaviour of the scaling functions. Thus, estimates of the inclusive top cross section at the LHC based on the threshold approximation alone are unreliable.

Results for bb production at HERA-B
In this section we present the inclusive cross section for bottom quark production at fixedtarget pp experiments, in particular, the HERA-B experiment. The energy of the proton beam at HERA-B is 920 GeV so that √ S = 41.6 GeV. Here the gg channel is dominant (about 70% of the total cross section), with the qq channel contributing the remainder. The gq and gq channels are again negligible, of the order of a few percent. Figure 13 shows that the η ∼ < 1 region is dominant in the convolution with the parton densities, Eq. (19).
In Fig. 24 we present our NLO and NNLL-NNLO results for the b quark production cross section in fixed-target pp interactions as a function of beam energy in the range 200-1200 GeV/c with m = 4.75 GeV. Comparing 1PI and PIM kinematics we find that, particularly at high energies, the NNLO predictions are different. Again, this can be attributed to increased sensitivity to the high-energy asymptotics of the scaling functions. In Table 3 we list the NLO and NNLO inclusive b quark production cross sections. We show the effects of both using the damping factor 1/ √ 1 + η and the approximate or exact (Eqs. (85)-(90)) scaling functions f (2,1) ij and f (2,2) ij on the NNLO results. We observe that the damping factor has a significant effect at small scales. The use of exact scaling functions also has a dramatic effect since it causes the NNLO bottom cross section in PIM kinematics to become negative at µ = m/2. The 1PI inclusive cross section at NNLO is significantly larger than the PIM cross section. Only the PIM NNLO results show a significant    reduction in scale dependence. Note also that for most options the scale uncertainty is larger than the kinematics ambiguity. In general, we see that the theoretical control over this cross section is still rather poor.

Conclusions
In this paper we have resummed threshold enhancements to double-differential heavy quark hadroproduction cross sections, both in single-particle inclusive and pair-invariant mass kinematics and in the qq and gg production channels. We used these resummed expressions to construct analytic approximations for the heavy quark cross sections up to NNLO, including the first three powers (NNLL) of the large logarithms. This involved including color-coherence effects and soft radiation attached to all the one-loop virtual corrections, the contribution of which we determined via matching conditions. We derived exact results for the scale-dependent terms at NNLO using renormalization group methods. We examined the magnitude of the NNLO correction to the inclusive cross section at the parton and hadron levels and studied its variations due to kinematics choice, factorization scheme choice, and changes in the renormalization/factorization scale for top quark production at the Tevatron and bottom quark production at HERA-B. We found that the scale dependence of the cross section after including the NNLO correction is typically substantially reduced, as expected on general grounds [59,60].
We now provide what are, in our judgement, reasonable values of the approximate NNLO top cross section at 1.8 and 2.0 TeV and the approximate NNLO bottom cross section at 41.6 GeV. These values, based on a combination of the results and their uncertainties, are obtained in the following procedure. For each energy, we take the AP+DF results for each kinematics choice, given in the tables in section 6. The scale uncertainty for 1PI or PIM is given by the maximum and minimum values of the cross sections in each row. The central value of the cross section is the average of the µ = m values in 1PI and PIM kinematics. Two errors are assigned to the resulting cross section. The first is due to the kinematics-induced ambiguity and is the difference between the central value and the µ = m values obtained in 1PI and PIM kinematics alone. The second is the weighted average of the scale uncertainties in the two kinematics, giving more weight to smaller scale uncertainties. Thus we obtain the following NNLO top production cross sections at the Tevatron, and σ tt (2.0 TeV) = 8.0 ± 0.6 ± 0.1 pb .
Recall that the first set of errors indicates the kinematics ambiguity while the second is an estimate of the scale uncertainty. Note that the scale uncertainty is considerably smaller than the kinematics uncertainty for top production. There are other sources of uncertainties in the total cross section such as parton distribution function uncertainties [61,62,63,64,65], ambiguities in the analytical expressions near threshold, as discussed in the previous section, and contributions from subleading logarithms [30]. Contributions from even higher orders may also be significant [30]. Due to the relative large uncertainties in the cross section we give one less significant figure in our combined results than in the tables of the previous section. We note that the central values for the cross sections do not change at all for top production and change only slightly for bottom production if we use any of the other options in the tables. The kinematics uncertainty also changes only slightly for the other options. The scale uncertainty is increased substantially, however, if we use the EX or EX+DF options. The combined results presented here reflect our subjective judgment -more unbiased results are presented in section 6.
In conclusion, we have provided, using soft gluon resummation techniques, improved finiteorder perturbative estimates for heavy quark production along with their uncertainties. We hope that our results will thereby allow more meaningful comparisons with measurements.

A NNLL matching terms for PIM kinematics
Here we collect the terms required for NLO matching in PIM kinematics. We first define the Born functions that occur in these terms. Recall that in our notation the superscript (n) indicates the power, (α s (µ R )/π) n , in the expansion of the corresponding quantity.
The qq channel 1PI Born function is The PIM Born function is where the subscript PIM indicates that the expressions for t 1 and u 1 in Eq. (14) should be used with β M replaced by β.
The gg channel 1PI Born function is The color factors C O , C K and, for later use, C QED , are defined as The PIM Born function is We again indicate that the expressions in Eq. (14) should be used for t 1 and u 1 with β M replaced by β. We also define for PIM kinematics Let us now turn to the one-loop S + MF cross sections, as defined in section 4. The MS scheme 1PI qq result is The MS scheme 1PI gg result is The MS scheme PIM gg result is B NNLL differential partonic heavy quark cross sections to NNLO In this appendix we collect all our results for the NNLL-NNLO differential partonic heavy quark cross sections. We first make some general observations. The NLO 1PI results agree with those in Refs. [32,33,54]. The NLO PIM results are consistent with Ref. [34]. The qq channel results contain terms that are antisymmetric under t 1 , u 1 interchange, a consequence of charge conjugation asymmetry of the initial state. The gg channel results are symmetric under t 1 ↔ u 1 . In the qq channel the Born function factors out as a whole for LL and NLL terms. In the gg channel the Born function function factors out as a whole only for LL terms. We also use the notation that the superscript (n) denotes the coefficient of (α s (µ R )/π) n in the expansion of the corresponding quantity.

B.1 The qq channel in 1PI kinematics
In the qq channel, Eq. (4), the Born cross section is with F Born qq (s, t 1 , u 1 ) given in Eq. (A.1). The MS one-loop NNLL corrections are qq 1PI (s, t 1 , u 1 ) .
The MS two-loop NLL corrections are To achieve NNLL accuracy one must add qq 1PI (s, t 1 , u 1 ) 2 ln(s 4 /m 2 ) s 4 where we suppressed the qq channel label on the soft anomalous dimension matrix elements. We note that at NNLL accuracy we have derived all NNLO soft plus virtual terms involving the scale, except for δ(s 4 ) terms involving single logarithms of the scale. We can go beyond NNLL accuracy and also derive these terms in the partonic cross section by requiring that the scale dependence in the hadronic cross section cancel out. These terms are qq 1PI (s, t 1 , u 1 )/F Born qq (s, t 1 , u 1 ). We have checked that these results are consistent with the exact expressions in Eq. (86) and with the expansion of the resummed cross section beyond NNLL accuracy as discussed in Ref. [30].
The DIS one-loop NNLL corrections are The DIS two-loop NLL corrections are To achieve NNLL accuracy for the DIS scheme at two loops we have to add to the NNLL MS terms

B.2 The gg channel in 1PI kinematics
In the gg channel, Eq. (5), the Born cross section is with F Born gg (s, t 1 , u 1 ) given in Eq. (A.3). The MS one-loop NNLL corrections are gg 1PI (s, t 1 , u 1 ) .
The MS two-loop NLL corrections are To achieve NNLL accuracy one must add where we suppressed the gg channel label on the soft anomalous dimension matrix elements. As we discussed in the previous subsection, we can also derive the NNLO δ(s 4 ) terms involving single logarithms of the scale. These terms are withT ′ (1) gg 1PI (s, t 1 , u 1 ) =T (1) gg 1PI (s, t 1 , u 1 )/F Born gg (s, t 1 , u 1 ). We have checked that these results are consistent with the exact expressions in Eq. (89) and with the expansion of the resummed cross section beyond NNLL accuracy.

B.3 The qq channel in PIM kinematics
In the qq channel, Eq. (12), the Born cross section is To achieve NNLL accuracy to two-loops one must add 12 Γ

+ 4 ReΓ
(1) As for 1PI kinematics, we note that at NNLL accuracy we have derived all NNLO soft plus virtual terms involving the scale, except for δ(1 − z) terms involving single logarithms of the scale. We can go beyond NNLL accuracy and derive these terms in the partonic cross section by requiring that the scale dependence in the hadronic cross section cancel out. These terms are The DIS two-loop NLL corrections are s d 2 σ To achieve NNLL accuracy for the DIS scheme at two loops we have to add to the NNLL MS terms s d 2 σ

B.4 The gg channel in PIM kinematics
In the gg channel, Eq. (13), the Born cross section is  To achieve NNLL accuracy to two-loops one must add 22 − C A + C A ln 2 √ ν 1 ν 2 s m 2 (15), keeping only terms that behave as ln(β). We checked that we find the same results for both kinematics. Our results are as follows.
For the qq channel in the MS scheme