Inclusive Heavy Flavor Hadroproduction in NLO QCD: the Exact Analytic Result

We present the first exact analytic result for all partonic channels contributing to the total cross section for the production of a pair of heavy flavors in hadronic collisions in NLO QCD. Our calculation is an essential step in the derivation of the top quark pair production cross section at NNLO in QCD, which is a cornerstone of the precision LHC program. Our results uncover the analytical structures behind observables with heavy flavors at higher orders. They also reveal surprising and non-trivial implications for kinematics close to partonic threshold.


Introduction
Thanks to the Tevatron, the production of heavy flavors at hadron colliders has nowadays become a very dynamic field of research. The imminent physics expected from the Large Hadron Collider (LHC) places even stronger emphasis on this class of processes. The very large production rates for both top and bottom at the LHC will alow studies of heavy flavors with precision so far seriously considered only for processes with light flavors. Accurate knowledge of top production is fundamental in Higgs and most Beyond the Standard Model as well as Standard Model studies [1][2][3][4][5].
Various sources of uncertainty affect the theoretical predictions for heavy flavor hadroproduction. The most important one is due to unknown higher order corrections; these are NNLO effects in the strong coupling or electroweak effects. Their reduction mandates higher order calculations.
Partial NNLO corrections to the total cross section have been presented in the literature [25][26][27]. They originate from a fixed order truncation of the all-order exponentiation of soft gluon logarithms [28][29][30][31][32]. We would like to emphasize, however, that such an approach to the fixed order cross section does not represent a controllable approximation in the sense of being derived from a consistent expansion in a small parameter. A priori, such partial corrections cannot be expected to approximate well the exact fixed-order result; an explicit example can be found in Fig.3 in Ref. [33] as well as in Section 4 below. A careful investigation of this point is beyond the scope of the present article. It will be provided elsewhere [34].
In a recent analysis, the authors of Ref. [35] point out that bound-state effects may have been underestimated in past studies of kinematics very close to threshold. They argue that this effect can have a significant impact on both near-threshold differential distributions and top mass measurements.
Another source of uncertainty in top production are initial state non-perturbative effects contained in the parton distribution functions. These effects have been discussed extensively in the literature; for a comparative study of the various pdf sets in top quark pair production see for example Ref. [26,32]. The importance of top quark pair production in constraining pdf sets and other LHC observables has been studied in Ref. [36]. Overall, pdf-related effects in the considered process are expected to be in the 3-4 percent range and are smaller than the missing higher order effects of interest to us in this work.
The present paper represents a needed step in the program [37][38][39][40][41][42][43] for the systematic derivation of the NNLO QCD corrections to the total heavy flavor production cross section. The phenomenological importance of the total pair-production cross section was established in the original pa-pers [8], where it was demonstrated that the effect of the NLO corrections is mainly on the overall normalizations and not on the shapes of differential distributions (however see Ref. [35] as mentioned above). It is natural to expect that this feature will persist at even higher orders, especially for distributions not subjected to too strong cuts. The reason for this expectation can be traced back to the non-trivial kinematics of the produced heavy quark pair starting from the leading order approximation.
Specifically, in this paper we derive the first exact analytic expression for the NLO QCD corrections to the total inclusive cross section for heavy quark pair production at hadron colliders. The results obtained are thus new, since previously only a numerically extracted fit was available. We comment on the quality of this numerical approximation in Sec. 4. Secondly, we aim at tuning a calculational approach that is capable of tackling the derivation of the NNLO corrections for this observable.
The calculation of the top quark pair production cross section at higher orders presents several significant challenges. First, the analytic structure of the total cross section is very complicatedand different from the one encountered so far in multiloop massless calculations. Indeed, radicals start to show up even at leading order. Our approach in addressing this problem is discussed in Section 3. Our findings (see sections. 4,5,6) not only confirm the expectations, but also demonstrate that a priori unforeseen analytic structures arise starting from NLO. Second, at NNLO one also has to tackle the problem of the large number of three-loop cut diagrams. Third, in the case of top quark pair production the value of the top mass is not small compared to the typical partonic center of mass energy. That in turn means that one has to obtain results, which are exact in the mass. In consequence, approaches based solely on an expansion in powers of the mass of the heavy flavor are not likely to be sufficient to cover the whole kinematic range. We illustrate this point by analyzing the newly derived exact NLO result in Sec. 4.
Finally, we would like to address the question of why we prefer an analytic approach and do not pursue a purely numerical one from the very beginning. One reason is that the NLO result, including terms subleading in the regularization parameter, is a prerequisite for the NNLO calculation. It is clearly desirable to have complete control over the lower order result. In fact, as we show in our following discussion, effects due to uncertainties inherent in numerical calculations can be sizable enough to require additional consideration. Our second motivation for pursuing analytical methods is based on the observation that the cross section we compute essentially represents a one-variable problem. Therefore, it seems reasonable to expect that at least most of the calculation should be feasible in terms of known functions and by using automatic methods. Our expectations have indeed turned out to be true. As a result we have gained deep insight into the analytic structure of cross sections as complex valued functions in the massive case.
The paper is organized as follows: in the next section we introduce some basic notations and definitions. In Sec. 3 we introduce our calculational method. Section 4 contains our main result, the analytic expressions for all partonic channels contributing at NLO. Sections 5 and 6, present detailed discussions of the analytic features of our results. Finally, in order to be as self-contained as possible, we collected various known formulae in the Appendix.

Setup
We calculate the NLO corrections to the total cross section for the process At partonic level, this involves the following reactions Following the established notations, we denote the cross section for the reaction Eq. (1) as where s = x 1 x 2 S is the partonic center of mass energy, S is the corresponding hadronic invariant and m is the pole mass of the heavy quark. The physical region is defined by s ≥ 4m 2 > 0. The scale µ denotes both the renormalization and factorization scales as appropriate. For simplicity, we do not distinguish between the two. The functions f i are the parton distribution functions for parton i in the (anti)proton.
The partonic cross sectionσ can be written as an expansion in the strong coupling. We will keep the traditional notation, where the cross sections are written in terms of the following dimensionless functionsσ where The partonic cross sectionσ contains three implicit scheme dependencies. First, it depends on the definition of the strong coupling α s . We work in terms of the standard MS coupling defined through where S ε = (4π) ε exp(−εγ E ) and β 0 = 11 6 C A − 2 3 T F n f is the first term of the QCD β-function (known by now up to the four-loop level [44,45]). The color factors in an SU(N)-gauge theory are C A = N, C F = (N 2 −1)/(2N) and T F = 1/2. Throughout this paper n f stands for the total number of flavors, which is the sum of n l light and n h = 1 heavy quarks. Switching between definitions of the coupling with n l + 1 (as adopted in this paper) and n l active flavors is done with the standard decoupling relations (also known up to the four-loop level [46,47], see Appendix A for the expression needed in our case).
Second, the partonic cross sectionσ depends on the definition of the mass m. As we mentioned above, we work in terms of the usual pole mass [48,49]. The third implicit dependence is on the scheme used to factor out the collinear (mass) singularities related to initial state collinear radiation. We adopt the standard MS factorization scheme defined througĥ where ρ = 4m 2 /s with 0 < ρ ≤ 1. The variable τ, ρ ≤ τ ≤ 1, stands for the corresponding (dummy) convolution variable. The functions σ kl are the "bare" (i.e. collinearly un-renormalized) hard scattering cross sections and the collinear counterterms Γ are given by where P (n,n l ) (τ) are the (n + 1)-loop (space-like) QCD splitting functions [50][51][52][53][54][55][56][57]. The superscript n l in P (n) and α s emphasizes that these quantities are evaluated in a scheme with n l -flavors.

Parton level evaluation
In this section we describe the evaluation of the partonic cross section σ i j (ρ, ε) appearing in Eq. (7). Following the method introduced in [58] we perform simultaneously all loop and real radiation integrations. With the help of the IBP identities [59,60] which we solve with the Laporta algorithm [61], the whole problem is mapped into a total of 37 master integrals.
The needed partonic cross sections are non-trivial functions of a single dimensionless variable ρ = 4m 2 /s 1 . As was established in the literature, the analytic structure of the result becomes rather complicated. In particular, even at LO radicals like β = √ 1 − ρ appear which suggests that the application of established analytic methods in this problem may be too cumbersome. To rectify the situation we perform the following change of variables This change of variables allows one to express the results for the partonic cross sections almost entirely in terms of harmonic polylogarithms (HPL) [62] of the variable x as well as rational functions of x and 1 ± x. Such analytic structures have been widely explored and can be dealt with in a straightforward manner. A particularly welcoming feature of this approach is that one can achieve a large degree of automatization in the solving of the differential equations for the masters with the help of existing software [63,64].
The reason why this change of variables is helpful can be seen as follows. The physical cross section for the production of two heavy particles in the final state has two physical singularities: the points m 2 = 0 representing the high energy (or massless) limit, and the kinematic threshold for heavy pair production s = 4m 2 . The third singularity, which is not in the physical region, is |m 2 /s| → ∞. The transformation (9) maps these three singularities into the points x = (0, 1, −1) which are precisely the singular points of the differential equations defining the set of HPL's.
As we will see in the following, however, the differential equations for a few of the master integrals posses singularities other than the three described above. These are the points s = m 2 , s = −m 2 , s = −4m 2 and s = −16m 2 . They are all outside the physical region 0 < 4m 2 /s ≤ 1. To the best of our knowledge, this is the first calculation of an observable with massive flavors where such pseudo-threshold-like singular points have been encountered. At NLO they appear only in the g + g → QQ + X subprocess starting from order O(ε 0 ). We discuss them in detail in Sec.5.
To completely fix the solutions for the master integrals one has to supply suitable boundary conditions. The behavior of the solution around any one of the singularities x = (0, ±1) can be used. In the cases where integrals form a system of two or three differential equations it turns out that one needs to supply information about the behavior of the solution at more than one singular point. It is also worth mentioning that the approach used in Ref. [65,66] to fix the boundary conditions by simply requiring regularity of the solution cannot be applied to the problem at hand since the masters are typically singular at all singular points. A detailed discussion of the most complicated master we encountered in this calculation is given in Sec. 6. Here, let us only note that while some of the boundaries required a direct calculation by Mellin-Barnes methods [67,68], where we used the MB package [69], many could be simply obtained by requiring the vanishing of the integral at threshold.
The evaluation of the convolutions in Eq. (7) also deserves a comment. A complication arises from the fact that while the dummy variable τ is the 'natural' variable for the splitting functions entering the collinear counterterms Γ, the same variable is not the natural one for the partonic cross sections σ. To write the latter in their natural variables one has to switch from the variable ρ to the conformal variable x defined in Eq. (9).

Results
At this point, we are ready to present our main result, namely the analytic expressions for the partonic cross sections at next-to-leading order. Our calculation has been performed with full account of the dependence on the number of colors N. Because of the substantial length of the formulae, however, we will present here only the result for N = 3. The complete results, valid for any SU (N) color gauge group are available in electronic form in [70].
4π f Figure 1: Error implied by the numerical integration of the cross section [6] in the quark-antiquark production channel against the exact result. The functions are defined in Eqns. (10,25). Figure 2: Error implied by the numerical integration of the cross section [6] in the quark-gluon production channel against the exact result. The functions are defined in Eqns. (11,25).
H(0, 0, 1, x) H(0, 1, 0, x) Most of the results above are expressed in terms of the standard harmonic polylogarithms  Figure 3: Error implied by the numerical integration of the cross section [6] in the gluon-gluon production channel against the exact result. The functions are defined in Eqns. (12,25). The plot on the right contains also six terms of the threshold (long-dashed line) and the high energy (shortdashed line) expansions.
(HPL) [62] of weight up to three and of the simple argument x, which implies that they can also be rewritten in terms of the usual polylogarithms of more complicated arguments but same weight.
As we already mentioned in Section 3 and further elaborate in Section 5, the presence of additional singularities beyond those naively expected makes it impossible to integrate the gluon-gluon cross section in terms of HPL's alone. Indeed, four additional functions are need. The first three of them are chosen as follows It is clear from these expressions that both F 1 and F 2 can also be expressed through the usual polylogarithms of up to weight three. Since the resulting expressions are very lengthy and cumbersome, yet obtainable by automated software, we prefer to stick to the integral representations which are easier to manipulate. Unlike F 1,2 , however, the presence of the square root in the integrand of F 3 makes it impossible to give a closed form result in terms of standard special functions in that case.
Besides being able to provide numerical values for the partonic cross sections at any given point, we would also like to perform threshold and high energy expansions. As far as HPL's are concerned this can easily be done with the help of [63]. As it turns out, F 1 , F 2 and F 3 can be expanded with little additional effort by simply integrating the expansions of the integrands term-by-term. The missing boundaries are trivial at threshold, since all integrals vanish there by their very definition. In the high energy limit, the asymptotic behavior necessary to complete this procedure is The origin of the last function, F 4 , is discussed at length in Section 6. We choose to write its integral representation in terms of the original mass variable ρ = 4m 2 /s, not x, (recall that the two are related through Eq. (9)) because this simplifies the arguments of the special functions. We have where The integrand is expressed through the complete elliptic function of the first, K, and second, E, kind defined as Obtaining the asymptotic behavior of F 4 in the two interesting limits is a substantially more complicated task than in the case of F 1 . . . F 3 . This is best done with the help of differential equations. Here we only reproduce the leading term of the expansions F 4 (β) = − 69 80 Having given all of the relevant expressions, we would now like to address the question of the actual precision of the fitting formulae derived by numerical integration in [6]. To that end, we plot the error defined as where f (1),NDE i j has been taken from [6], as function of the velocity, β, over the full range from the threshold β = 0 to the high energy limit β = 1. As can be seen in Figs. (1,2,3) the error exceeds the promised 1% only in the vicinity of a zero of the cross section correction. This is of course expected, since there is no particular need or phenomenological interest to find the exact position of the zeros only to have an improved fitting formula, and the authors of [6] have not used this information to obtain their approximations. Overall, as can be concluded from these plots, the quality of the fits in the case of f (1) qq and f (1) gq is excellent. However, the error almost approaches the 1% bound in the case of f (1) gg . While some loss of precision is expected given how complicated a function the gluon-gluon cross section is, we will see in the following that this slightly larger uncertainty is of relevance to the threshold behavior of the cross section.
The analytic results presented above make it possible to obtain exact expansions in both the massless and the threshold limits. As far as practical applications at colliders are concerned, the threshold behavior is more important. We give it here with full account of the color factors and T F set to its customary value of 1/2 gg can be found in Appendix A. The singular term proportional to 1/β (Coulomb enhancement) as well as the coefficients of the logarithms log(8β 2 ) (soft-gluon enhancement) have been known already from [6].
The β-independent constant terms are important for the quantitative understanding of heavypair production at threshold. Numerical approximations were first derived in [6]. Recently the au-thors of Ref. [35] independently extracted these terms from earlier calculations [71,72] on quarkonium production at hadron colliders. The comparison of the results extracted in [35] with the one derived from the calculation of [6] showed that the two results were incompatible with each other, both for qq and gg scattering. The numerical size of the discrepancy was found to be significantly larger than the numerical uncertainty inherent in the results of Ref. [6].
With the help of our calculation this discrepancy is now resolved and the conclusions are quite intriguing. In the qq channel the discrepancy turned out to be due to a missing contribution in the quarkonium calculation, see Ref. [35] for details. The situation is quite different in the case of the gluon-gluon channel. In this reaction our calculation agrees with the analytic result extracted in Ref. [35]. Thus, the source of the discrepancy is solely due to an accidental magnification of the numerical uncertainty of the calculation in Ref. [6]. Indeed, the numerical counterpart of the β-independent constant term in Eq. (27) is according to [6] 768π 7 a gg 0 ≃ 37.25 .
The difference between the two results is about 7%, even though the fitting formula for the cross section from Ref. [6] is more accurate than 1%, see Fig. 3. This is not entirely surprising as the following arguments show. Close to threshold the term proportional to this coefficient in f (1) gg is damped by a factor of β, whereas the dominant contribution comes from logarithms of the velocity and Coulomb enhanced constant terms. Thus its error is less relevant for the final precision. For larger values of β other terms dominate and their errors conspire to make the result sufficiently precise. An analysis showing that this mechanism indeed works, i.e. changing the value of a gg 0 by up to 10% leaves the total cross section within its 1% bounds, has been performed in [35].
This difference represents a powerful example for the need of analytic results in certain cases. It will be very interesting to further study how this effect propagates in the soft-gluon resummed cross section. Work is in progress [34].
Let us finally turn our attention to the high energy limit of the cross sections. It is well known that the t-channel gluon exchange leads to a constant behavior of the quark-gluon and gluon-gluon cross sections in this regime. We derive The constant in Eq. (30) is in perfect agreement with the numerical value obtained in [6]. The relation between the two cross sections as given in the second equation above has been explained there as well.
The quality of deeper expansions around β = 0 and β = 1 is shown in Fig. 3. With six terms in both cases, most of the range is covered satisfactorily. In fact, due to large numerical cancellations in the vicinity of β = 0, it is better to use the expansion rather than the exact formula. Let us stress that these cancellations are a problem only in the case of numerical evaluation. Higher precision in the numerics always allows to retrieve the correct result.

Singularities of the differential equations
As already mentioned in Sec.3, besides s = 4m 2 , m 2 = 0 and s = 0, we have observed four additional singular points of the differential equations of the master integrals. These additional singularities characterize the integrals depicted in Fig. 4, which contribute exclusively to the gluon-gluon cross section.
The presence of additional singularities in the differential equations may have two different implications. One possibility is that, in principle, a different choice of master integrals may remove these singularities if they do not have some deeper (physical) meaning. However, if these singularities persist in the amplitude, it immediately follows that it would not be possible to express the final result through harmonic polylogarithms of the x variable alone.
We study these alternatives in the following. The main idea behind our analysis is that, although we work with cut graphs, the homogenous parts of the corresponding differential equations are the same as for ordinary Feynman integrals, which means that if the corresponding Feynman integral has a singularity then the cut graph may have it as well.
It turns out that the singularity at s = −16m 2 occurring in the graph of Fig. 4 a) is the one easiest to explain. Indeed, a quick look at the Landau equations shows that this is the leading singularity of the Feynman integral, when all lines go on shell. As such, it can be classified as an anomalous pseudo-threshold, because it occurs in the Euclidean region. How to integrate the differential equations in this case is explained in Section 6.
The second singular point, which is also rather easy to explain, is s = −4m 2 . Let us take a closer look at Fig. 4 b). The graph is crossed in the forward direction because of the exchanged outgoing momenta. Its existence is clearly due to the possibility to switch the identical external gluons. If we treat the graph as a normal box graph with arbitrary external momenta, then it would have a physical singularity at t = 4m 2 . The forward limit of the crossed graph corresponds here to the Mandelstam u = 0, and not t = 0. This in turn implies that t = −s and the usual branch point modeled by a logarithm is mapped to −s = 4m 2 . Interestingly, this argument further implies that this singularity can occur in any crossed box graph (in the sense explained above), because the heavy quark line always has to form a closed loop, and therefore the singular point t = 4m 2 is always present. Interestingly, the seemingly simplest graphs in our problem have the singularity that is hardest to explain. The singularity of Fig. 4 c) at s = m 2 does not correspond to any leading singularity in the sense of Landau equations. If we kept the external momenta as free, then there would be three master integrals for this topology and mass distribution. Their differential equations do not have more singularities than those of the one-loop graph without the massless internal line. It is only when we reach the forward direction that the three integrals collapse to one and generate the aforementioned singularity in the process. Fortunately, Fig. 4 c) and d) are linked. In fact going from one to the other while forgetting the cut, can be done by the change s ↔ t. Since in the forward direction t = −s, the singularity at s = m 2 of c) implies a singularity of d) at s = −m 2 .
A closer examination of the final result for the gg cross section Eq. (12) shows that all four singular points described above are present in the final result in some form or another. Indeed, Figs. 4 a), b), c) and d) correspond to F 4 , F 3 , F 1 and F 2 respectively. Whereas s = −16m 2 for F 4 and s = −m 2 for F 2 are usual branch points of logarithmic type, s = −4m 2 is only a square root type branch point of F 3 , and s = m 2 is not a branch point of F 1 at all. This is exactly the opposite of what would happen if we had to do with an amplitude, because there, all points in the Euclidean region would have to be regular.
As a final remark we would like to stress that while at next-to-leading order these additional singularities are only present in the gluon-gluon channel, we anticipate that they will be present also in the quark-anti-quark channel at the next-to-next-to leading order.

The most difficult master integrals
Of all 37 master integrals occurring in our problem, two stand out as particularly difficult. They both correspond to the topology and mass distribution depicted in Fig. 5 and contribute to the gluon-gluon cross section exclusively. As usual in such cases, integration-by-parts (IBP) identities make it possible to choose any powers of the propagators and irreducible numerators, as long as one of the integrals cannot be reduced to the other (independence with respect to IBP's). We make the following choice where D i are Feynman propagator denominators of the lines which are not cut, p 1 and r are defined in Fig 5, whereas k and l are the momenta of the remaining cut lines also shown in Fig. 5. The "+" sign of the δ-functions defines the energy flow of the cut lines, as usual.
One can show that both M 1 and M 2 are finite (with or without the cut). Moreover, they enter the cross section with coefficients regular in the dimensional regularization parameter ε, which means that we only need to know their value at ε = 0. In terms of the dimensionless variable ρ = 4m 2 /s, they satisfy the following system of differential equations Although their exact form is irrelevant for our forthcoming discussion we reproduce here the nonhomogenous terms H 1 and H 2 Let us rewrite the system as one second order equation for M 1 . We obtain where The homogenous part of Eq. (38) defines a Gauss hypergeometric function. A change of variable to k = √ −4ρ, further simplifies the equation The two solutions of the homogenous part are where K is the complete elliptic function of the first kind defined in Eq. (21).
As it stands, the solution Eq. (42) requires analytic continuation in the physical range, because k is actually imaginary. Using the relation we choose to set At this point, we only need to compute the Wronskian of the solutions in order to be able to present a complete solution to our problem. Using the fact that the derivative of the complete elliptic function of the first kind is where E is the complete elliptic function of the second kind defined in Eq. (22), we obtain the following Wronskian In the above equation we have also used the Legendre relation Finally, the solution to Eq. (38) is The choice of the upper bound above allows to easily find the integration constants C 1 and C 2 .
After obtaining M 2 from Eq. (34), and using the fact that both integrals have to vanish at threshold, it turns out that C 1,2 = 0.
Before we turn to some general conclusions based on the result we just obtained, let us stress that in this case the right choice of integration variable has proven to be invaluable. Indeed, had we tried to solve the differential equations in our conformal variable x, it would have been impossible to identify the differential equations as hypergeometric.
A quick inspection of Eq. (48) shows that it is impossible to express this integral with polylogarithms or a generalization thereof (defined as nested integrals over some elementary functions). This statement is in particular true of Goncharov's polylogarithms, and is a simple consequence of the definition of the elliptic integrals. In general, the presence of functions with such properties in master integrals does not preclude the possibility for their cancellation in the final cross section. However, as the results in Section 4 demonstrate, these non-polylogarithmic functions persist also in the final result.
The fact that the total cross section requires a substantially extended set of functions is probably not that surprising, because of the never trivial phase space integration. After all, the two masters considered in this section correspond to a three particle cut. However, the analysis we have just performed implies that the same functions will be necessary to solve the master integrals for the fermionic contributions to the two-loop amplitude for heavy flavor production in gluon fusion. Examples of Feynman diagrams that contain the integral from Fig. 5 without the cut are given in Fig. 6. In this argument it is crucial that the only place where the presence of a cut plays a role are the non-homogenous terms H 1 and H 2 (and therefore the function G). In consequence, the result for the two-loop integral will be given by the very same Eq. (48). This means that it will be substantially more difficult to obtain an analytic result in this case than it was in the quark annihilation channel discussed in [42]. Finally, let us stress that even though elliptic integrals have been encountered previously in Quantum Field Theory (see for instance a recent study in [73]), we believe that our observation in this particular problem is an important lesson on complications associated with massive QCD processes. Whereas in the massless case ordinary polylogarithms are sufficient, with masses present on internal lines one will encounter elliptic functions and possibly other yet unexpected objects.

Summary
In this work we have calculated for the first time the exact analytic expressions for the next-toleading QCD corrections to all inclusive cross sections for heavy quark pair production at hadron colliders. The most important phenomenological application is top quark pair production at the Tevatron and LHC.
The same corrections have been computed in approximate form by several groups using more traditional methods based on numerical phase space integration. Our calculation utilizes completely independent methods where both real and virtual integrations are performed simultaneously and in analytic form. We find complete agreement with the earlier calculations; the comparison to our exact result shows that the numerical uncertainties in the qq and gq reactions are indeed very small. While much larger when compared to the case of the other two reactions, the uncertainty in the gg reaction does not exceed the promised accuracy of 1%.
The motivation for the present calculation was manifold. First, it is a step in the ambitious program for the calculation of the NNLO corrections to top quark pair production at the LHC. Indeed, UV and collinear renormalization require the computation of terms subleading in the corresponding regulator, or in the case at hand, beyond d = 4. Second, one would like to derive and test methods capable of tackling the very demanding NNLO calculation. Third, since the resulting complexity is driven by the underlying analytic structures, the current calculation allows one to gauge the true level of the complexity of the NNLO calculation. We find for the first time that new, a priori unexpected analytic structures do appear even at NLO. What is more important, we are able to pinpoint their specific origin. Indeed, at NLO we have to compute a total of 37 master integrals and only seven of them have an analytic form differing from what is known from purely massless calculations. Therefore, the conclusion that we draw from our present work is that the analytic methods applied here alow a large degree of automatization and also seem suitable for the calculation of (at least) the bulk of the NNLO result.
Finally, the analytic results derived in this paper have important implications to the NLO/NLL (and beyond) phenomenology of top quark pair production at the Tevatron and LHC, especially for the production close to the partonic threshold. In particular, we demonstrate that although the uncertainty in the previous numerical evaluations of the gluon-gluon cross section is within 1% of the exact result, the extraction of the so-called "constant" terms that are very important at threshold leads to a much larger numerical difference, one that is in fact as large as 7%.
Our findings for the threshold behavior of the cross section have to be carefully explored, especially in light of the fact that the progress and improvements in the top quark pair production cross section in the last ten years or so were based on refinements of the behavior of the cross section at threshold. A detailed investigation of this effect will be presented elsewhere [34].  26) and (27), and setting T F to its customary value of 1/2