Top-quark pair production near threshold at LHC

The next-to-leading order analysis for the cross section for hadroproduction of top quark pairs close to threshold is presented. Within the framework of non-relativistic QCD a significant enhancement compared to fixed order perturbation theory is observed which originates from the characteristic remnant of the 1S peak below production threshold of top quark pairs. The analysis includes all color singlet and color octet configurations of top quark pairs in S-wave states and, for the dominant configurations, it employs all-order soft gluon resummation for the hard parton cross section. Numerical results for the Large Hadron Collider at $\sqrt{s} = 14 $TeV and $\sqrt{s} = 10 $TeV and also for the Tevatron are presented. The possibility of a top quark mass measurement from the invariant mass distribution of top quark pairs is discussed.


Introduction
At the CERN Large Hadron Collider (LHC) the major part of top quarks are produced in pairs. Due to the experience gained at the Fermilab Tevatron [1] and the huge amount of top quarks to be produced at LHC the reconstruction of top quarks with good accuracy will be possible [2,3]. A significant fraction of top quark pairs will be produced close to threshold. Thus a dedicated analysis of the production cross section in this region is required which is best performed within the framework of non-relativistic QCD (NRQCD) [4,5].
The production of top anti-top quark pairs close to the kinematical threshold has received much attention in the context of precision measurement of top quark properties at a future International Linear Collider (ILC). Theoretical calculations and dedicated experimental analyses have demonstrated that a precise extraction of the top quark mass, its width and the strong coupling constant is possible [6,7] at the ILC. The complete next-to-next-to-leading order (NNLO) predictions are available since many years [8]. (For earlier work see e.g [9][10][11][12].) Partial next-to-next-to-leading logarithmic (NNLL) [13,14] and next-to-next-to-next-to-leading order (NNNLO) [15][16][17] predictions were evaluated more recently.
In contrast to the linear collider, where the physical observable is the total cross section as a function of energy, at the hadron collider one considers the invariant mass distribution of the top quark pairs. Since the expected uncertainty is significant larger than the one anticipated at a linear collider a next-to-leading order (NLO) analysis is probably sufficient. The calculation of the cross section within the NRQCD framework contains as building blocks the hard production cross section for a top quark pair at threshold and the non-relativistic Green's function governing the dynamics of the would-be boundstate. Both ingredients are available in the literature since many years. In particular, the hard cross section for threshold tt production can be found in Refs. [18,19]. In Ref. [18] the NLO formulae were derived for quark or gluon initial states and a quarkonium in a J PC = 0 −+ color singlet state, plus possibly a parton. The general case, with the heavy quark system (QQ) in S -wave singlet/triplet spin state, and color singlet/octet configuration is given in Ref. [19], together with the corresponding results for P-waves. The results of Refs. [18,19] were presented for stable boundstates. For unstable wide resonances it is convenient to describe the bound state dynamics through a Green's function.
velocity. Finally, we perform a soft gluon resummation which enhances the cross section by a few per cent.
Our paper is organized as follows: In the next Section details of the formalism used for the calculation of the NLO cross section are provided. The effects of initial-state radiation and the hard contribution are discussed in Section 3 and the soft gluon resummation is performed in Section 4. The properties of the Green's function are summarized in Section 5. In Section 6 the building blocks are combined and numerical results for the invariant mass distribution are presented. Theory uncertainties due to scale variation and unknown higher order corrections are estimated. Summary and conclusions are presented in Section 7.

The production cross section
Let us denote the (quasi) boundstate of a top and anti-top quark with spin S and angular momentum L by T ≡ 2S +1 L [1,8] J where the superscripts [1] and [8] denote the singlet and octet color states. The production rate is obtained from the production cross section of a top quark pair with invariant mass M 2 ≡ (p t + pt) 2 and its evolution to a quasi boundstate described by the non-relativistic QCD. The former is a hard QCD process at a distance ∼ 1/m t and thus computable within the conventional perturbative expansion in α s .
The long-distance effects responsible for the formation of a narrow boundstate are described by the squared wave function at the origin |Ψ(0)| 2 or, in the language of NRQCD, by the matrix elements Here N s = 2S +1 and N c = 1, 8 denote the number of spin and color degrees of freedom, respectively. We are interested in the differential distribution dσ/dM which, for narrow resonances with mass M n , is proportional to δ(M − M n ). For wide resonances, the case under consideration, it is convenient to convert the factor describing the sum over individual resonances into the non-relativistic Green's function 1 with G(M + iΓ t ) ≡ G [1,8] ( r = 0; M + iΓ t ) being the Green's function at zero distance for the nonrelativistic Schrödinger equation discussed below. Since the typical momentum scale governing the nonrelativistic top quark system m t (with m t 2 ≡ M + iΓ t − 2m t , and being the velocity of top and anti-top quarks) is in the perturbative regime, and the large top quark width Γ t introduces an additional cutoff scale √ m t Γ t , the Green's function can be evaluated perturbatively. As stated above the present paper is concerned with the production of top quark pairs near threshold, thus restricted to states with L = 0, i.e. T = 2S +1 S [1,8] J . The contributions to the invariant mass distribution with higher angular momentum are at least suppressed by 2 , and thus of higher order (beyond NLO).
In order to obtain experimentally measurable quantities at a hadron collider the partonic differential cross section dσ i j→T /dM is convoluted with the luminosity function where i, j refer to partons inside the hadrons P 1 and P 2 with the distribution functions f i/P 1 and f j/P 2 . The dependence on the factorization scale µ f cancels in combination with the one contained in dσ i j→T /dM. The differential cross section can thus be written as As usualŝ and S denote the partonic and the hadronic center-of-mass energy squared, respectively, and τ =ŝ/S . The lower limit of the τ integration is given by ρ = M 2 /S . The partonic differential cross section dσ i j→T /dM consists of a factor F that is evaluated in perturbative QCD, and can be deduced from Refs. [18,19], and a second factor, the imaginary part of the Green's function G [1,8] M ImG [1,8] where the superscript of the Green's function refers to the color state of T . Eqs. (4) and (5) constitute our master formulae, which contain several scales and various physics contributions of different origin in factorized form. In particular, the soft dynamics of the parton distribution and real radiation is contained in the convolution of F i j→T with the parton luminosity, the boundstate effects are described by G. Note that at NLO the Green's function G [1,8] (M + iΓ t ) and the convolution of F with the parton luminosity (L ⊗ F) are individually independent of the renormalization scale µ r . Thus we can discuss the two parts separately in the following two Sections. Furthermore, it is simpler to assess the uncertainties for the individual contributions.
Let us at this point make a comment concerning the validity of Eq. (5), which makes use of the NRQCD expansion assuming ≪ 1, thus being limited to the threshold region. For larger invariant masses conventional perturbation theory is applicable (see Refs. [22][23][24] and Refs. [25][26][27][28] for recent compilations of the total cross section and [29] for a proposal to measure the top-quark mass from the shape of dσ/dM). In the transition region the predictions from both methods are expected to coincide, as will be discussed below (c.f. Fig. 4).
The coefficients C h originate from the hard corrections to the production process. The functions A c contain the real corrections with collinear parton splitting from one of the initial partons i, j, and are governed by the Altarelli-Parisi splitting functions, A nc originates from non-collinear real emission. These individual contributions are manifest already in Ref. [18] and the appendix of Ref. [19], and will be listed in the following. Note, that in Eq. (6) we have split off the factor (1 + (α s /π) C h ), which we attribute to hard corrections and thus treat as a multiplicative factor to the terms in square brackets.
In Tab. 1 we collect all processes of the type i j → T X at NLO which contribute in our analysis and list the corresponding normalization factors N i j→T . Note that the production of a spin triplet color singlet state 3 S [1] 1 via gq or qq scattering is zero up to and including NLO. This is because in these channels the heavy quarks are produced through gluon splitting g * → tt, which is only possible if the tt is in an octet state.
The coefficients C h are non-vanishing only for the processes which are present also in lowest order [18,19]: where β 0 = (11/3)C A − (4/3) n f T F and C F = 4/3, C A = 3, T F = 1/2, n f = 5. The last term in C h [qq → 3 S [8] 1 ], arising from non-decoupling of the top quark in the gluon propagator, has been observed and discussed in Ref. [20], see also footnote 3 on page 73 of Ref. [20]. For the other processes hard corrections are of higher order, thus C h is zero at NLO: The function A c is conveniently expressed using Altarelli-Parisi splitting functions P i j (z) introduced below [18,19] where the conventional plus-distribution 2 was employed to regularize the singularity at z = 1. The splitting functions P i j (z) are given by 2 The plus-distribution follows the prescription 1 0 dz ln n (1−z) is an arbitrary test function which is regular at z = 1. It is related to the ρ-prescription used in Ref. [19] by The functions A nc are obtained from the non-collinear contributions. For spin singlet states we have where ] is singular at z = 1, and regularized by the plus-prescription. For spin triplet states one obtains The function A nc [qq → 3 S [8] 1 ] is also defined with the plus-prescription. The leading singular behavior of A nc is given by A nc (z) [8] 0 and qq → 3 S [8] 1 . In the soft limit its behavior is insensitive to the details of the boundstate and only depends on its color configuration.
It is instructive to discuss the relation between the normalizations of different processes leading to the same boundstate. For instance, the normalization N i j→T (see Tab. 1) for the process gq → 1 S [1,8] 0 X is fixed by gg → 1 S [1,8] 0 , because in the collinear limit this cross section factorizes into the corresponding LO process and the P gq splitting function. As a consequence the cancellation of the factorization scale dependence happens among the gg and gq initiated reactions. Similarly, the normalization of gq → 3 S [8] 1 is fixed by qq → 3 S [8] 1 . In contrast, the processes qq → 1 S [1,8] 0 and gg → 3 S [1,8] 1 are forbidden at LO, hence the corrections have to be collinearly finite. In comparison to Ref. [20] the combinations A c + A nc include terms that vanish in the limit z → 1. Furthermore subprocesses that appear for the first time in O(α 3 s ) were neglected in Ref. [20]. The relative size of these terms will be adressed below.
Let us now start the numerical analysis. The partonic cross sections have to be convoluted with the parton distribution functions (PDFs) in order to arrive at the hadronic cross section. We use the CTEQ6.5 [30] set for the PDFs and take α (5) s (M Z ) = 0.118, m t = 172.4 GeV and √ S = 14 TeV as input values. The running of α (5) s (µ r ), which is the input for the partonic cross sections, is evaluated with the help of RunDec [31], using the four-loop approximation of the β function. This leads to α (5) s (µ r ) = (0.1077, 0.09832, 0.09050) for µ r = (m t , 2m t , 4m t ). Furthermore we identify renormalization and factorization scales (µ f = µ r ).
As stated above, the cross section factors into the convolution L ⊗ F and the Green's function. To discuss the relative importance of the various contributions individually the results for the subprocesses without the factor ImG(M + iΓ t )/m 2 t are given in Tab. 2. Note that color-singlet tt production is dominated by by gg → 1 S [1] 0 . Color-octet production is dominated by gg → 1 S [8] 0 plus a 25% contribution from qq → 3 S [8] 1 . The size of the remaining subprocesses (neglected in Ref. [20]) amounts to five to ten percent and is strongly scale dependent. The variation of µ (recall µ = µ f = µ r ) between m t and 4m t leads to changes of L ⊗ F by ±3% and ±7% for the total singlet and octet production, respectively. In these channels the real radiation of partons contains large logarithmic contributions in the NLO corrections. In combination with the rapidly varying parton luminosity these logarithms make up for a major part of the numbers quoted in Tab. 2. The origin of these large logarithms can be traced to the singular behavior of the cross section near z ≈ 1, regularized by plus-distributions. There exists well established technology for the resummation of these large logarithms to all orders in perturbation theory. We will address this issue next.

Soft gluon resummation
The parton channels, which exhibit enhancement due to soft gluon emission are gg → 1 S [1] 0 , gg → 1 S [8] 0 , and qq → 3 S [8] 1 (see Eqs. (9) and (12)). The relevant logarithms are contained both in A c (from initial state radiation) and A nc (from FSR) and read for the three leading processes: where To investigate the size of the threshold logarithms, we evaluate the contribution of the factorized hard scattering contribution convoluted with the PDFs, i.e. L ⊗ F separately for the three contributions: tree-level, singular and regular terms. (The hard corrections (1 + (α s /π) C) are common to all). The threshold enhanced contributions are defined in Eq. (13) and correspond exactly to the terms included in Ref. [20], while regular terms correspond to the remainder of A c + A nc in Eqs. (9) and (12) The three lines correspond to µ = µ f = µ r = (m t , 2m t , 4m t ). We note that in all three cases the contribution of the threshold enhanced terms from Eq. (13) is large, although the regular terms in the case of qq → 3 S [8] 1 are of the same order. Technically the matching applied in Ref. [20] corresponds to neglect all terms which vanish exactly at threshold that is for z = 1, i.e. Eqs. (11) and (12) of Section 3. The regular terms in Eq. (14), which have not been accounted for in the recent analysis of Ref. [20], are of the same order as the NLO sub-processes as given in Tab. 2.
Threshold resummation proceeds conveniently in Mellin-space. To that end we calculate the Mellin moments with respect to z = M 2 /ŝ according to Then, the Mellin-space expression for the threshold enhanced terms listed in Eq. (13) read (see also [32,33]) where we have kept all dominant terms in the large-N limit and neglected power suppressed terms of order 1/N. γ E is the Euler-Mascheroni constant (γ E = 0.577215 . . .).
The resummed expressions (defined in the MS-scheme) for the individual color structures of the hard cross sections F of Eq. (6) are given by a single exponential in Mellin-space (see e.g. Refs. [33][34][35]) where F (0),N i j→T denotes the tree level term in Eq. (6) and the exponents are commonly expressed as ln ∆ N i j→T = ln N · g 1 i j (λ) + g 2 i j→T (λ) + . . . , where λ = β 0 α s ln N/(4π). To next-to-leading logarithmic (NLL) accuracy the (universal) functions g 1 i j as well as the functions g 2 i j→T are relevant in Eq. (18), see Ref. [25] for the extension to NNLL accuracy. Explicit expressions are where the full dependence on µ r and µ f has been kept. The gluonic expressions g 1 gg and g 2 gg→T are obtained with the obvious replacement A (i) q → A (i) g . The perturbative expansions of the anomalous dimensions are universal and well-known. We have [36] and all gluonic quantities are given by multiplying A (i) q by C A /C F . We also give explicit results for the matching functions g 0 i j→T in Eq. (17), For phenomenological applications [37,38] of soft-gluon resummation at the parton level one introduces an improved (resummed) hard cross section F res , which is obtained by an inverse Mellin transformation as follows, The NLO results can also be found in Tab. 2.
Here F NLO i j→T is the standard fixed order cross section at NLO in QCD, while F N i j→T | NLO is the perturbative truncation at the same order in α s obtained by employing Eq. (16). That is to say that for the matching we have fully expanded all formulae consistently to O(α s ). This adds the hard coefficients C h of Eq. (7) to the results Eqs. (16) and (21). In this way, the right-hand side of Eq. (22) reproduces the fixed order results and resums soft-gluon effects beyond NLO to NLL accuracy.
In Section 6 we employ Eq. (22) for phenomenological predictions by performing the inverse Mellin transform numerically. To that end, one should note that the treatment of the precise numerical matching to the exact NLO hard cross section is a matter of choice since different schemes lead only to differences which are formally of higher order. We have found that the application of the resummed result is well justified when the kinetic energy of the top-quark pair is a few GeV or less, see e.g. Ref. [25], where the precise numerical value is not important. Another issue concerns the constant terms in Eq. (21) which are sometimes modified to include formally sub-leading (but numerically not insignificant) terms, see for instance Ref. [37,38]. As just explained, in the present analysis we adopt the minimal approach, i.e. we apply Eq. (21) (including the hard coefficients C h of Eq. (7)) and account for all regular terms in Eq. (14) through matching to NLO.
In Tab. 3 we compare the fixed-order NLO and resumed result of the convolution L ⊗ F. One observes an enhancement up to about 10% depending on the process.

Boundstate corrections
Let us next discuss the boundstate corrections. As mentioned above, the convolution of F i j→T with the parton luminosities provides the normalization of the differential cross section, while its shape is mainly determined by the non-relativistic Green's function. The latter describes the long-distance evolution of the top quark pair produced near threshold. The kinematics of the produced top quark pair is nonrelativistic, and the dynamics is governed by exchange of potential gluons leading to the formation of quasi-boundstates. The corresponding potential is given at NLO by with C [1] = C F = 4/3 and C [8] = C F − C A /2 = −1/6, and a 1 = (31/9)C A − (20/9) T F n f .
The color-singlet Green's function feels an attractive force, the color-octet Green's function is governed by repulsion and thus does not develop a boundstate. They are both defined as the solutions of the Schrödinger equations For the Green's function at zero-distance, the NLO result is known in a compact form [39] (see also [14]) G [1,8] (M + iΓ t ) ≡ G [1,8] ( r = 0; M + iΓ t ) = C [1,8] with κ ≡ iC [1,8] Here  25) is an expansion around the LO boundstate poles and thus has multiple poles of a form G ∼ |Ψ (0) n (0)| 2 /(M (0) n − M) k (k = 1, 2 at the NLO). However, resummation of this multiple poles into single poles is straightforward and well-known. We refer to Ref. [39] for further details.
In Fig. 1 we show the imaginary parts of the color singlet and color octet Green's functions in the threshold region. As input we use m PS t = 170.1 GeV, which to NLO accuracy corresponds to m t = 172.4 GeV [1], and Γ t = 1.36 GeV [40][41][42]. At NLO the Green's function is separately renormalization scale invariant and we are free to chose µ r independent from the hard process. A well-motivated physical scale is µ s = m t C F α s (µ s ) = 32.21 GeV which corresponds to twice the inverse Bohr radius. The corresponding α s value used in Fig. 1 is α (n f =5) s (µ s ) = 0.1401. It has been observed that the color-singlet Coulomb Green's function has a well-convergent perturbative series for this scale choice [43].
In order to see the effect of Coulomb resummation, we plot for both color states three lines: the full Green's function (solid line) and the expansion of G in fixed order up to O(α s ) with and without top quark width (dashed/dotted). The upper three lines in Fig. 1 correspond to the color singlet case and the lower three to the color octet one. The color-singlet Green's function shows a pronounced peak which corresponds to the tt resonance below 2m t , while for color octet there is no enhancement. Note that the curve for the full octet Green's function is very close to the one-loop expansion (taking into account the finite top quark width). Thus for the color octet state the Coulomb resummation effect is negligible. In addition, one more line (dash-dotted) for the color-singlet Green's function is plotted including the NNLO Coulomb potential, which is useful to estimate yet unknown boundstate corrections to the NLO colorsinglet Green's function. As input value we again adopt the PS top quark mass [44] given above. Note that in the absence of full NNLO result for the Green's function and hard correction, this improved Green's function would not be sufficient for a full NNLO prediction. Nevertheless, the difference between solid and dash-dotted curves gives an indication of the intrinsic uncertainties of the Green's function, which is roughly 10%. i + α s C [1,8] In the zero-width limit (iΓ t → +i 0), the color-singlet curve for the expansion exhibits a step of height α s C F /8 (for M → 2m t ), and the color-octet curve formally becomes negative for ≤ −α s C [8] π/2 which corresponds to M − 2m t < 0.23 GeV. Both for the singlet and octet case the fixed order result without Γ t the imaginary part of the Green's function vanishes below 2m t . The qualitative difference between the solid and the short-dashed curves will be reflected in the comparison of our final results for the invariant mass distribution with the prediction based on a fixed order calculation: for the color-singlet curve we observe a sizable excess in the region below the nominal threshold up to roughly 5 GeV above. In the color-octet case, as a consequence of the relative smallness of C [8] , the prediction follows roughly the Born approximation. Although the color-octet Green's function is significantly smaller than the singlet one, the relatively large hard scattering factor L ⊗ F for 1 S [8] 0 plus 3 S [8] 1 , which exceeds the one for the singlet case by roughly a factor four, quickly over-compensates the effect of the Green's functions.
In the present paper we use the analytical result of the Green's function, which includes the α s correction (i.e. the second term in the square brackets of Eq. (23)) by means of the Rayleigh-Schrödinger perturbation approach. In Ref. [20] a numerical solution to Eq. (24) has been employed, which resums the α s corrections to all order. The numerical solution is more stable against scale variation and applicable over a wide range of µ r . However, the difference between the two approaches is below 2% and formally of higher order. Extensive studies on higher order effect to the color singlet Green's function exist in the literature (see, e.g., Refs. [8,43]), including different implementations of the Green's function. From the experience collected in the linear collider studies on tt production, we expect rather large corrections from the variation of µ r for the color singlet Green's function of about 20% which is significantly bigger than the estimate from the NNLO Green's function mentioned above. In contrast to the color-singlet case the higher order corrections to the color octet Green's function are expected to be unimportant since there is no resonance enhancement and the color coefficient C [8] is small.

Invariant mass distribution
We are now in the position to combine the results of the preceeding Sections and discuss the cross section for the invariant top quark distribution.
In Fig. 2 the invariant mass distributions for LHC ( √ S = 14 TeV) is shown for the three dominant processes. The bands reflect the scale variation of the convolution L ⊗ F which for the color singlet case amounts to roughly ±1%. The reduction as compared to Tab. 2 and Fig. 2 is due to a compensation of the µ dependence after including the sub-leading NLO processes. Note, however, that the corresponding Green's function shows an uncertainty due to the renormalization scale variation of about 20% which is well-known from top quark production studies gg → 1 S 0 Figure 2: Invariant mass distributions for leading subprocesses: gg → 1 S [1,8] 0 (blue and light green, respectively) and qq → 3 S [8] 1 (green). For each process the bands take into account scale variation of the hard cross sections.
in e + e − collisions, consistent with the difference between solid and dash-dotted curves in Fig. 1 and thus not discussed in Fig. 2. This pattern is also evident from Fig. 3, where all production channels as listed in Tab. 2 are included. The width of the bands is obtained from varying renormalization and factorization scales in the hard cross section as described above. The additional uncertainty from the Green's function, which we estimate 20% for the singlet and below 5% for the octet case, is not included.
As expected, for M < 2m t the production of tt pairs is dominated by the singlet contribution. However, for M > 2m t one observes a strong raise of the octet contributions, in particular of gluon induced subprocess which for M ∼ > 2m t +5 GeV becomes even larger than the corresponding singlet contribution. For the color-octet case the scale dependence of the hard scattering amounts to ±7%. Considering the threshold behavior as shown in Figs. 2 and 3 it is clear, that the location of the threshold is entirely governed by the behavior of the color singlet (S -wave) contribution. Thus, as a matter of principle, determining the location of this step experimentally would allow for a top quark mass measurement, which is conceptually very different from the one based on the reconstruction of a (colored) single quark in the decay chain t → Wb. In fact, much of the detailed investigations of tt threshold production at a linear collider were performed for this particular relations between the location of the color singlet quasi-boundstate  pole of tt and the top quark MS-mass. The absolute normalization of the cross section is also sensitive towards electroweak corrections [45][46][47][48][49] which are of the order of 5% close to threshold. For example, the difference between corrections from a light (M h = 120 GeV) and a heavy (M h = 1000 GeV) Higgs boson amounts to roughly 6% [48].
In Fig. 4 the prediction for dσ/dM based on NRQCD is compared with the one obtained from a fixed order NLO calculation for stable top quarks which is obtained using the program HVQMNR [50]. As expected from the comparison of solid and dotted curves in Fig. 1, the two predictions overlap for invariant masses around 355 GeV. Above 355 GeV relativistic corrections start to become important. From this comparison we find an additional contribution to the total cross section for tt production of roughly 10 pb, which could become of relevance for precision measurements. Note that the band of the NRQCD-based prediction only contains the uncertainty from the scale variation of L ⊗ F whereas the one of the Green's function (which can reach up to 20%, see Section 5) is not shown.
The analysis of this work has concentrated on the threshold region and is applicable for M up 360 GeV at most. However, it is obvious, that the overall shape of dσ/dM will be distorted and the mean M shifted to smaller values, which might affect the global fit of dσ/dM. In Fig. 5 we present for comparison the NLO prediction for dσ/dM in the wide range up to 700 GeV.  The distribution reaches quickly its maximum of 3.3 pb/GeV at around 390 GeV and then falls off slowly. It is remarkable that its value at 370 GeV is already not too far from the maximum of the curve and the threshold modifications thus affect a sizeable part of the distribution.
Although the most detailed top quark studies will be performed at the LHC at an energy of 14 TeV, a sample of top quarks has been collected at the Tevatron in proton anti-proton collisions at 1.96 TeV. Furthermore the first LHC data set will be taken at 10 TeV. For this reason we give the results for these two cases, in Figs. 6 and 7. The cross section in Fig. 6 has the same characteristic shape as the one in Fig. 3, however, the absolute size is considerably smaller. As expected, the enhancement at threshold is significantly less pronounced for Tevatron where the colour singlet contribution is very small.
Our analysis confirms the findings of Ref. [20], however, the numerical results for the cross sections as presented in Fig. 3 are slightly higher than the corresponding corrections of Ref. [20] which is due to the combined effect of the soft-gluon resummation, the inclusion of the NLO sub-processes and the different matching to full QCD.

Summary
A NLO analysis of top quark production near threshold at hadron colliders has been performed. The large width of the top quark in combination with the large contribution from gluon fusion into a (loose bound) color singlet tt system leads to a sizable cross section for masses of the tt system significantly below the nominal threshold. A precise measurement of the M tt distribution in this region which is dominated by the color singlet configuration could lead to a top-quark mass determination which does not involve the systematic uncertainties inherent in the determination of the mass of a single (colour triplet) quark. Furthermore, also the shape of the differential distribution dσ/dM is distorted and the mean M shifted towards smaller values.
The effects of initial state radiation as well as boundstate corrections are taken into account in consistent manner at NLO. As compared to Ref. [20] we include the completeŝ dependence in the matching condition and also implement all NLO sub-processes. We observe a partial numerical cancellation between these two effects leading to similar predictions as Ref. [20]. Furthermore we perform a soft-gluon resummation and thus include the dominant logarithmically enhanced higher order terms. This last step stabilizes the prediction. However, it enhances the cross section at most by 10%.  The effects are more pronounced at the LHC with top production being dominated by gluon fusion and less relevant in proton-antiproton collisions with top quarks dominantly in color octet states. Considering the threshold region (say up to M tt = 350 GeV) seperately, an integrated cross section of 15 pb is obtained, which should be compared to 5 pb as derived from the NLO predictions using a stable top quark and neglegting the binding correction. Within this relatively narrow region the enhancement amounts to roughly a factor three and a significant shift of the threshold. Compared to the total cross section for tt production of about 840 pb (obtained using fixed-order NLO accuracy for µ = m t , see, e.g., Ref. [25]), the increase is relatively small, about 1%. However, in view of the anticipated experimental precision of better than 10% these effects should not be ignored.  At the Tevatron qq → 1 S [8] 0 dominates the cross section, and luminosity for gg channels is small, thus the boundstate peak is buried by color-octet production.

Note added
While this article was finished an analytic evaluation of the total cross section at NLO accuracy appeared [51], which has been used in Ref. [20] to clarify the existence of a non-decoupling top quark effect overlooked in Ref. [19] (see footnote 3 on page 73 in Ref. [20]).