Virtual amplitudes and threshold behaviour of hadronic top-quark pair-production cross sections

We present the two-loop virtual amplitudes for the production of a top-quark pair in gluon fusion. The evaluation method is based on a numerical solution of differential equations for master integrals in function of the quark velocity and scattering angle starting from a boundary at high-energy. The results are given for the renormalized infrared finite remainders on a large grid and have recently been used in the calculation of the total cross sections at the next-to-next-to-leading order. For convenience, we also give the known results for the quark annihilation case on the same grid. Outside of the kinematical range covered by the grid, we provide threshold and high-energy expansions. From expansions of the two-loop virtual amplitudes, we determine the threshold behavior of the total cross sections at next-to-next-to-leading order for the quark annihilation and gluon fusion channels including previously unknown constant terms. In our analysis of the quark annihilation channel, we uncover the presence of a velocity enhanced logarithm of Coulombic origin, which was missed in a previous study.

to provide a boundary condition in the form of high precision values of the integrals at a single point. Inspired by Refs. [24,25], a point in the high-energy range has been chosen for this purpose. The asymptotics of the integrals in this limit have been derived using Mellin-Barnes techniques. The second problem is related to singularities of the differential equations, which cause substantial problems. In the case of the quark-annihilation channel amplitudes the use of higher numerical precision was sufficient to provide numerical values within some acceptable kinematical domain. Unfortunately, the gluon fusion channel is substantially more demanding both in the determination of the asymptotics of the integrals, and in the treatment of numerical instabilities. The techniques that we have applied in this case are documented in this work.
The presentation of numerical results is always a non-trivial issue. If the functional dependence of the amplitudes is smooth enough, one might consider giving a grid of values. This is indeed what we have done. One must, however, remember that due to well-known singularities (mostly Coulomb and collinear), the amplitudes diverge in some corners of the phase space. A result in the form of a grid is not practicable there. Instead, we provide expansions beyond the borders of the grid, both close to the threshold and far away from it. Finally, even though the quark-annihilation amplitudes are known since Ref. [12], we reproduce them here in the same format as those of the gluon-fusion channel.
The two-loop amplitudes presented here have already been used for the calculation of total cross sections. Moreover, they will be used in the near future for the evaluation of differential distributions, and will certainly serve as a benchmark for on-going analytic calculations. There is yet another application of our results, which concerns the threshold expansion of the total cross sections. In Ref. [26], we have discussed how to derive such an expansion from soft-gluon factorization using virtual amplitudes and soft functions. In this work, we will apply these methods and derive the leading threshold effects including quarkvelocity independent terms. This is an extension of the results of Ref. [27]. Interestingly, while the study presented in the latter reference aimed at the derivation of all the terms singular in the velocity of the top-quark, our new result shows a discrepancy in the form of a logarithm of the velocity, which is due to Coulombic effects. We will elucidate the origin of the omission.
The paper is organized as follows. In the next section, we introduce our notation, and discuss ultraviolet and infrared renormalization of the amplitudes. Subsequently, we describe the evaluation of the high-energy asymptotics of the master integrals and the numerical/semi-numerical solution of differential equations for arbitrary kinematics including expansions at singular points. We continue with the presentation of the results for the amplitudes including benchmark numerical values, and threshold expansion. Finally, we present the threshold expansions of total cross sections. The main text is closed with conclusions and followed by appendices containing renormalization constants and anomalous dimensions, as well as cross section expansions for arbitrary color. The numerical results for the amplitudes on a phase space grid, and their high energy expansions are attached to the electronic submission of this paper.

Preliminaries
We consider the virtual amplitudes for two processes: the gluon-fusion channel heavy-quark pair-production g(p 1 ) + g(p 2 ) → Q(p 3 , m) +Q(p 4 , m) , (2.1) and the quark-annihilation channel heavy-quark pair-production q(p 1 ) +q(p 2 ) → Q(p 3 , m) +Q(p 4 , m) . (2.2) In both cases, the final state heavy-quarks (top-quarks) are on-shell, p 2 3 = p 2 4 = m 2 . The results presented in the next sections for the first process are new, whereas those for the second process are an improvement, as far as the kinematical range is concerned, over those given in Ref. [12]. The amplitudes can be expressed in terms of modified Mandelstam invariants and θ is the scattering angle. The ratio of the mass and the center-of-mass energy is expressed in terms of the velocity of the top-quark This is a natural variable for current hadron collider applications, where it is known that the bulk of the top-quark pair-production cross section comes from the range of moderate velocities. The high-energy limit is condensed to β ≈ 1, while the production threshold is to be found at β ≈ 0. The bare on-shell amplitudes admit a perturbative expansion, of which only the first three terms are relevant (2.6) The subscript g or q specifies the initial state, while the ket-notation signifies a vector in color and spin space. We have also indicated the dependence on the parameter of dimensional regularization with d = 4 − 2 space-time dimensions. The external degrees of freedom are treated with Conventional Dimensional Regularization (CDR), in particular the gluons have d − 2 degrees of freedom. As usual, we are interested in color and spin summed squared amplitudes. The two-loop contributions are then given as a product of the Born and two-loop amplitudes, and may be written as an expansion in the number of colors, N c , of a general SU(N c ) group, and the number of closed light-, n l , and heavyquark, n h , loops. The heavy quarks in the loops are assumed to have the same mass as the external top-quarks. We write The factor 2 difference between the prefactors in the equations has no deeper meaning, but is taken over from [24,25]. The amplitudes are renormalized according to where Z g , Z q and Z Q are the on-shell wave-function renormalization constants for gluons, and light-and heavy-quarks respectively. A prefactor has been introduced in order to keep the amplitudes dimensionless and avoid unnecessary γ E − ln(4π) terms. The bare mass is just m 0 = Z m m, while the bare coupling constant is which corresponds to the MS scheme with n f = n l + n h active flavors, if the loop integrals are calculated with the measure d d k/(2π) d . The necessary renormalization constants can be found in Appendix A. Notice that our results are always given with µ = m. Renormalization with n l + n h active flavors has been assumed in every previous calculation of two-loop virtual amplitudes for heavy-quark pair-production. However, once we are interested in a hadronic cross section, we have to reassess the question of the number of active flavors. Indeed a natural scale for top-quark pair-production is set by the topquark mass. The bulk of the cross section comes from the regime, where non-logarithmic top-quark mass effects are not negligible, since the top-quarks themselves are not ultrarelativistic. It is thus more appropriate to calculate cross sections with n l active flavors and not resum the mass logarithms in either the coupling constant or the parton distribution functions. Amplitudes with n l active flavors are obtained by decoupling in α s , which amounts to the substitution α where the decoupling constant, ζ αs , has a perturbative expansion in the strong coupling constant and can be found in the Appendix A.
Although the ultraviolet renormalized amplitudes still contain divergences of infrared origin, the structure of these divergences is completely understood at the two-loop level [28][29][30][31][32][33][34]. We, therefore, define the finite remainder of the amplitudes through infrared renormalization s , m, µ, , (2.12) where we have underlined that the amplitude on the right hand side of the equation has been obtained by decoupling of the heavy quark from the running of α s . The MS renormalization constant Z Mg,q is a matrix in color space and has a non-trivial dependence on the kinematics {p} = {p 1 , ..., p 4 }, and by the same on the masses {m} = {0, 0, m, m} of the external partons. It can be derived from the differential equation where we do not specify the initial state anymore, and the color space matrix anomalous dimension is given by up to two-loops by [33] (2.14) The summations run over massless (indices i, j, k) and massive (indices I, J, K) partons, with the notation (i, j, ...) denoting unordered tuples of different indices. The color operators T a i act on the color indices of the respective partons. If the particle is a gluon carrying a color index c, we have (T a ) bc = −i f abc , assuming the result has been projected on color index b. Similarly, for an outgoing quark (or incoming anti-quark) the generator is (T a ) bc = T a bc , whereas for an incoming quark (or outgoing anti-quark) the generator is (T a ) bc = −T a cb . The kinematic dependence is contained in s ij = 2σ ij p i · p j + i0 + , where the sign factor σ ij = +1 if the momenta p i and p j are both incoming or outgoing, and σ ij = −1 otherwise. For massive partons there is p 2 I = m 2 I , v I = p I /m I , and cosh β IJ = −s IJ /2m I m J . It was noted in Refs. [32,33] that the triple-color-correlators in the third and fourth line of the anomalous dimension do not contribute to top quark pair production amplitudes' divergences. The general case was analyzed in Ref. [26] with the conclusion that these terms never contribute to color and spin summed amplitudes, as long as the latter are real. The anomalous dimensions' coefficients relevant to the present work are given in Appendix A.
With the two definitions of ultraviolet and infrared (finite remainder) renormalized amplitudes, we are ready to proceed with the description of the computational methods and results.

Methods
The two-loop amplitudes for heavy-quark pair-production are expressed through 726 Feynman diagrams in the gluon-fusion channel, and 190 diagrams in the quark annihilation channel. Due to the structure of the QCD vertices, the topologies present in the quarkannihilation case are a subset of those present in the gluon case. Using the Laporta algorithm, the occurring integrals are reduced to a set of 422 masters, out of which only 145 are needed for the quark-annihilation case. Based on experience, we consider integrals with less than six propagators as easy. This leaves 212 difficult six and seven line integrals to evaluate. In our work, we do not use any external input, i.e. we do not rely on integrals calculated by others. Nevertheless, 38 of the difficult integrals have been evaluated by one of us in Ref. [12]. The remaining number is further reduced by the fact that 89 integrals can be obtained from others by a t ↔ u transformation. Thus, the final number of new integrals evaluated in this work is 100. Thanks to the work Ref. [25], 17 of these were at least known in the high-energy limit.
As explained in the introduction, we chose to work as proposed in Ref. [25] and solve differential equations for the master integrals with numerical methods starting from a boundary at high-energy. Since we are dealing with four-point amplitudes with a single mass scale, the integrals, once stripped of their global mass dimension by appropriate rescaling, depend on two dimensionless variables. This means that the functional dependence of the integrals is fully specified by two systems of homogeneous linear first order partial differential equations which can be obtained by taking derivatives in the parameters and reducing the resulting integrals with integration-by-parts identities to the original masters. The elements of the Jacobi matrices, J (m 2 ) and J (t) , are rational functions of m 2 /s, t/s and . We require a solution for the master integrals in the form of Laurent expansions in . The original equations are, therefore, transformed into by expanding all quantities in to sufficient order. The tilde denotes the coefficients of the respective expansions. The solution of the system Eq. (3.2) is obtained by choosing a path, possibly complex, in the parameter space and solving the resulting ordinary differential equation in z. This can be done either numerically or as a power-logarithmic series in z. In practice, we have proceeded as follows 1. We have determined analytically the first few terms of the high-energy expansion of the master integrals. The results are a power-logarithmic series in m 2 /s, with coefficients, which are exact in t/s. In order to obtain our results, we have used Mellin-Barnes techniques, and in particular relied heavily on the MB package [35].
In a few cases, we recovered the exact dependence on t starting from the limiting behavior at t = 0 and using differential equations in t. In order to obtain the boundary condition, we performed a double expansion of the Mellin-Barnes representation of a given Feynman integral in m 2 and t. The resulting Mellin-Barnes integrals, which were pure numbers were evaluated with very high precision and resummed with the PSLQ algorithm [36]. In simpler cases, we have used the XSummer package [37] for resummation.
2. In a next step, we have obtained deep high-energy expansions using the differential equations in m 2 and the boundary conditions from the previous step. These expansions were used to derive high-energy expansions of the amplitudes, and to obtain high precision numerical values of the master integrals at small mass.
3. Using the numerical values determined in the previous step, we have solved the differential equations in m 2 and t along contours in the complex plane. To obtain the solution, we have used the software from Ref. [38] with improvements to handle higher precision numbers [39].
4. Around β = 0, we have obtained, with the help of the differential equations, deep power-logarithmic expansions in β for fixed values of the scattering angle. These expansions were generated from unknown boundary coefficients, which were determined by matching the expansion to the numerical solution from the previous step at a point, at which both the expansion and the numerical solution provide high precision. This method can be used at any singular point, and allows to avoid numerical instabilities of the differential equations.

Results
The results obtained with the methods of the previous subsection fall into three kinematical domains: threshold, "bulk", and high-energy. The "bulk" domain covers moderate β values and is given purely numerically on a large grid. The sampling values of β are chosen as in [1][2][3][4], i.e.
and β 80 = 0.999. This covers the range of values available at LHC @ 8 TeV. The dependence on the scattering angle is described through where the x i correspond to the 21 sampling points of the Gauss-Kronrod integration rule of order 10, and can be obtained with any major algebraic/numeric computer system, e.g. Mathematica. The Gauss-Kronrod rule is an efficient deterministic rule for smooth functions, which also provides an error estimate by sampling every second point of the rule with appropriate weights. This specific choice of the cos θ points has been made, because a first aim was to provide very precise values of the contributions of the amplitudes to the total cross sections. As explained in the previous section, our final results are given for the finite remainders of the renormalized virtual amplitudes with n l active flavors. The results are normalized in such a way that they directly give the contribution to the respective cross section. In particular we define In both channels, gluon-fusion and quark-annihilation, there iŝ whereσ (2) (β) is the two-loop contribution to the cross section. Since our results might also be used for other processes, such a b-quark pair-production, we keep the dependence on the number of active flavors exact, and further decompose The values of A (g) and A (q) on the grid are attached to this publication in electronic form. We give five significant digits for each phase space point. This is sufficient for any practical applications, even if interpolation errors and possible cancellations with the oneloop-squared contributions are taken into account. In order to illustrate the complexity of the amplitudes, we plot them in Fig. 1.
As expected, the functions are very smooth. In the case of gluons, we observe additional enhancements at β ≈ 1 and cos θ ≈ ±1. For low scattering angle, these are due to diagrams of the general form depicted in Fig. 2, which contain a t-channel heavy-quark propagator. At tree-level, for instance, we have The behavior at large angles is obtained by t ↔ u symmetry. As in previous publications on the subject, we also provide high-precision values for the color coefficients, Eqs. (2.7,2.8), of the amplitudes at a benchmark point. They can be found in Tab. 1 for the renormalized two-loop virtual amplitudes with n l + n h active flavors. The numbers for the quark-annihilation channel have been presented previously  in [12]. The table also shows, which coefficients are already known analytically. In Tab. 2, we present the results at the same phase space point, but for the finite remainders with n l active flavors. Values of the color coefficients defined in Eqs. (2.7,2.8) of the renormalized two-loop virtual amplitudes with n l + n h active flavors for the gluon-fusion and quark-annihilation channels at m 2 /s = 0.2, t/s = 0.45 with µ = m. The entries with explicit 0.0 vanish because of the scale choice. The gray shaded values are only known numerically, whereas the remaining ones have been obtained analytically in [19][20][21][22][23]33].
While the grid covers a substantial part of the phase space, it cannot completely cover the high-energy and threshold domains, since the amplitudes are singular in these limits. For these two cases, we provide expansions. The leading behavior of the high-energy expansions for the renormalized amplitudes has been first determined in Ref. [24,25]. We have calculated the first three terms of the expansions of the bare two-loop amplitudes analytically and converted them into finite remainders. The results are attached to this   The question of the validity of the high-energy expansions at small/large scattering angle deserves a more careful study. We discuss the small scattering angle case only, since the opposite limit is analogous. First, we note that at cos θ = 1, there is t/s = 1/2(1 − β) > m 2 /s, with t/s → m 2 /s for m 2 /s → 0. From this, we conclude that the ratio m 2 /t should be considered of the order of unity in the high-energy low-scattering-angle limit. The convergence of the expansions depends, therefore, on the coefficients of the series. By inspection, we note that the expansions seem to have a non-zero radius of convergence at cos θ = 1.
Independently of these considerations, in the high-energy small/large scattering-angle case, cross sections are dominated by real radiation, which means that the issue of the behavior of the virtual amplitudes in this region is moot.
In the threshold region, the amplitudes are dominated by singularities coming from potential interactions between the heavy quarks. Beyond the grid, sufficient precision is obtained already with the leading behavior up to terms of O(β 0 ). For the application of our results to threshold expansions of cross sections, we are also interested in decomposing the gluon-fusion channel amplitudes into color structures. In general, there are three SU(3) invariant color structures linking two incoming gluons with color indices a and b, with a final state heavy-quark with color index c and anti-quark with color index d: 1) color singlet, 1, given by δ ab δ cd ; 2) symmetric color octet, 8 S , given by d abe T e cd , with d abc = 2 Tr(T a T b T c + T c T b T a ); and 3) anti-symmetric color octet, 8 A , given by if abe T e cd . For each of these structures, denoted by α, we introduce a color projector, P α , such that P 2 α = P α , and P 1 + P 8 S + P 8 A = 1. The expansions of the finite remainders of the renormalized two-loop amplitudes with n l active flavors at µ = m read  +C F −8 + 6 ln 2 − 6 ln 2 2 + 7π 2 12 − 8 9 n h + n l − 10 9 + 2 ln 2 where the Born matrix elements near threshold are given by 14) The formulae for the expansions of the dominant channels contain constants, which we have only determined numerically. They are given in Tabs. 3,4. Notice that several coefficients in the second table (of fermionic and leading-color bosonic origin) could be determined analytically using the results of Ref. [19,20]. This requires, however, to add analytic results for the infrared counter-terms that transform renormalized amplitudes into finite remainders. In the gluon channel, analytic results are only known for the sum of the color-structure projected amplitudes. It is thus not possible to obtain any coefficients presented in Tab. 3 in analytic form with the currently available results from the literature.

Threshold expansions of cross sections
Our results for the two-loop virtual amplitudes can be used to obtain the leading threshold behavior of partonic top-quark pair-production cross sections at next-to-next-to-leading order. When β ≈ 0, cross sections are dominated by terms of the form β × 1/β n ln m β with n > 0 and/or m > 0. While the first factor of β is a phase space suppression, the enhancements by positive powers of ln β and 1/β are due to emissions of soft gluons and non-relativistic potential interactions between the quark and the anti-quark. At nextto-next-to-leading order, the coefficients of these singular terms have been determined in Ref. [27]. We would like to extend this analysis to include terms with n = m = 0, which are velocity independent with respect to the Born cross section. The threshold expansion including both singular and constant terms can be obtained from factorization as explained in Ref. [26]. According to the latter publication, close to threshold, a cross section for a given initial state can be written asσ H α are called hard functions, and are obtained by expanding in β the partonic cross sections obtained exclusively with the finite remainder of the virtual amplitudes projected onto the color configuration α. Therefore, H α do not contain any real-radiation effects. S α are called soft functions, and are given by cross sections for emission of gluons and light-quark pairs from eikonal lines representing the external partons of the hard process in the color configuration α. The convolution is performed in the energy of the soft radiation. The color configurations, α, must correspond to irreducible representations of the SU(3) group, in order for this simple form to be valid. The cross section expansions are thus labeled by the color configuration. Accordingly, we introduce the perturbative expansionŝ where we only consider the initial states ij = qq, gg, and the possible color configurations are α ∈ {1, 8} for the quark-annihilation channel, while α ∈ {1, 8 S , 8 A } for the gluon-fusion channel. In this section, we work exclusively with n l active flavors, which implies that s . We will again neglect the scale dependence and set µ = m in the subsequent expressions.
Using the matrix elements Eq. (3.14), we obtain the following threshold expansions for the dominant color-projected Born cross sections in the gluon-fusion channel As was pointed out in [40], it is a coincidence that the anti-symmetric octet cross section is sub-leading in the threshold region. This fact cannot be justified with symmetry arguments (the often quoted Landau-Yang theorem does not apply to gluons). There iŝ Finally, the quark channel receives octet contributions onlŷ Notice that while singlet contributions do not vanish beyond leading order in this channel, they are not enhanced at threshold and are of O(β 3 ). We are interested in the threshold expansions of the next-to-next-to-leading order contributions,σ (2) ij,α , which take the form The coefficients of the singular terms of the β-expansions, c (n,m) ij,α with n > 0 and/or m > 0, are in principle known [27], although a contribution has been omitted in the analysis of c (0,1) qq , as we explain in the following. For future reference, we define C (2) ij,α = c for all but the ij = gg, α = 8 A case. For the latter, we further define where the right hand side has been expanded to O(β 0 ). The coefficient C The purpose of this section is to obtain C (2) ij,α . However, we will first determine a missing contribution to c (0,1) qq . To this end, we return to the derivation of the coefficients of the β-singular terms in Ref. [27], which proceeded with an additional factorization of potential effects. These can, in fact, be classified into two categories. One category corresponds to the case, when a quark-anti-quark system is created, and subsequently strongly interacts through effective potentials, dominated by the Coulomb potential. An example diagram can be found in Fig. 3 a). This category gives a complete description of the singlet case, and was the only one taken into account in Ref. [27]. However, in the case of an octet final state configuration there is another category of enhanced contributions. Indeed, the pair can then also annihilate, as shown in Fig. 3 b). Here, the Coulomb exchange between the virtual heavy quarks leads to a logarithmic singularity as we will demonstrate below. Unfortunately, since this type of enhancements has not been taken into account in Ref. [27], the cross section expansions given there do not contain all the singular terms. The additional logarithm contributes to the quark-annihilation channel only. The reason for the absence of enhancement in the gluon-fusion channel is that the diagram of Fig. 3  b) can only occur in the anti-symmetric octet configuration, where the s-channel gluon is produced from a triple-gluon vertex. At next-to-next-to-leading order, this diagram is contracted with the Born amplitude, which is sub-leading in the anti-symmetric octet, and the contribution to the total cross section is only of O(β 3 ln β).
In principle, we could determine the correction to the singular expansion terms directly from Eq. (4.1) after substituting the expansions of the virtual amplitudes from the previous section. However, the missing term due to Fig. 3 b) can be obtained with simple arguments based on unitarity and analyticity, which helps to clarify the origin of the logarithmic enhancement. We start with the definition of the vacuum polarization for gluons i Π ab µν (q) = q , (4.9) and its decomposition by transversality Π ab µν (q) = (g µν q 2 − q µ q ν ) δ ab Π(q 2 ) . (4.10) With these definitions, it is easy to see that = Π(s) × . (4.11) In fact, the vacuum polarization insertions could be Dyson-resummed leading to the replacement . (4.12) By unitarity, the production cross section in the octet configuration is proportional to Im Π(s), which means that the result containing annihilation contributions to all orders is given by Im 1/(1 − Π(s)). Returning, however, to our fixed-oder analysis at next-to-nextto-leading order, we note that we only have to consider a single diagram, Fig. 3 b). The vacuum polarization contribution in this diagram will be denoted by Π (2) C (s). As long as we are only interested in Coulomb exchanges, the vacuum polarization close to threshold is expressed through the non-relativistic Coulomb Green function [41] with appropriate color-factor changes to take into account the octet configuration. Nevertheless, we will need much less information here. Indeed, if we consider the imaginary part of Π where, in the first approximation, we have used the fact that the imaginary part of the vacuum polarization is related to the cross section, which receives a Coulomb correction in the form of the Sommerfeld factor, πα/β/(1 − exp(−πα/β)), expanded to first non-trivial order, where now α = (C F − C A /2)α s . The second approximation follows from the wellknown (textbook) expression for the one-loop vacuum polarization function Π (1) (s). The β suppression of the imaginary part due to the vanishing phase space volume at threshold is canceled by the Sommerfeld 1/β enhancement. The expression thus tends to a nonvanishing constant at threshold. This is the reason for which the next-to-leading order cross section for top-quark pair-production also tends to a non-vanishing constant. At this point, we have only determined the imaginary part of the vacuum polarization in a particular limit, which does not allow us to use dispersion relations to obtain the real part. However, it is still possible to use unitarity and analyticity for the same purpose. Indeed, the singularity at β = 0 is a branch point modeled as usual by a logarithm. An appropriate variable to express the occurrence of the cut is the non-relativistic energy of the system, E = √ s − 2m ≈ mβ 2 , where the last approximation is valid for β ≈ 0. The imaginary part across the cut is obtained from (4.14) Thus, the fact that the imaginary part tends to a constant together with the presence of a cut in the non-relativistic energy starting from E = 0, allows us to write the final result for the contribution to the cross section (4.15) Let us now apply the factorization formula Eq. (4.1). For the case at hand, we writê whereσ (1), V, f in ij,α andσ (2), V V, f in ij,α are cross sections evaluated with the finite remainder of the virtual amplitudes at the next-to-and next-to-next-to-leading orders. Notice that the latter require the one-loop-squared corrections [42][43][44], besides the two-loop amplitudes of the present publication. In order to obtain the β-expansion of the one-loop-squared finite remainders, we have used our own results from [4]. The soft functions, S (i) ij,α , on the other hand, can be found in Ref. [45] for the singlet, and in Ref. [26] for the octet configurations of the final state. They are only affected by the choice of the initial state through the Casimir invariants of the respective representations of the gauge group.
Results forσ ij, α with exact dependence on the number of colors and active flavors can be found in Appendix B. Here, we reproduce the results with N c = 3 and n l = 5. For the gluon-fusion channel, there is [27] σ ( where the constant is obtained by combining the contributions from the three color configurations where the last number here and below corresponds to n l = 5. Combining these numbers leads to C (2) gg = 503.664 − 29.9249 n l + 0.142857 n l 2 = 357.611 . The correct value of the coefficient at n l = 0 differs from that from Ref. [4] by about 50%, which fits within the error estimate from the latter publication. A partial cancellation of this coefficient by the term proportional to n l is responsible for the even larger difference of the values at n l = 5.
In the quark-annihilation channel we have This formula differs from the one presented in Ref. [27], in the coefficient in front of ln β, by the amount given in Eq. (4.15). In the case of top-quarks, i.e. with n l = 5, the numerical effect of the correction is tiny. It amounts to a 2.5% reduction of the coefficient, which was quoted to be 528.557 in Ref. [27]. This is important, as the results for the exact total cross sections in Ref. [1] have been parameterized with respect to the threshold expansion. It turns out that the difference in the coefficient of ln β is invisible within the numerical precision of the results of Ref. [1]. For the constant term, we obtain The two results are compatible at the 10% level, which is exactly the uncertainty quoted in Ref. [1].

Conclusions and outlook
With this publication, we have completed the numerical analysis of two-loop amplitudes for heavy-quark pair production in the quark-annihilation and gluon-fusion channels. We have demonstrated that numerical methods based on differential equations lead to high precision results in the relevant kinematical range. We have also provided expansions in the threshold and high energy limits, which can be used to obtain reliable numbers for any phase space points. These results should also be viewed as benchmarks for future analytic evaluations. In any case, they will constitute the basis for the calculation of differential distributions for top-quark pair production, which is the subject of current work.
Besides amplitudes, we were also able to provide threshold expansions of partonic cross sections improving over previous studies. Our results can be incorporated into resummed predictions for total cross sections. Once virtual amplitudes become known completely analytically, it will be possible to give fully analytic results for the velocity independent terms of the threshold expansions. We have provided formulae, which should make such an exercise straightforward.