Two-Loop Leading Color Corrections to Heavy-Quark Pair Production in the Gluon Fusion Channel

We evaluate the two-loop QCD diagrams contributing to the leading color coefficient of the heavy-quark pair production cross section in the gluon fusion channel. We obtain an analytic expression, which is valid for any value of the Mandelstam invariants s and t and of the heavy-quark mass m. Our findings agree with previous analytic results in the small-mass limit and with recent results for the coefficients of the IR poles.


Introduction
The top quark is the heaviest elementary particle known to date. Its large mass of 173 GeV makes it a crucial tool for the study of the electroweak symmetry breaking mechanism and for the search of signals of physics beyond the Standard Model. The study of the properties of the top quark at the Tevatron during the last 15 years provided precise measurements [1] of the top-quark mass (with a relative error of 0.75%) and of the total top-quark pair production cross section (with a relative error of about 10%). Owing to the by-now large Tevatron dataset, measurements of differential distributions in top pair production are becoming possible [2,3]. A new set of precise experimental measurements is expected to be obtained by the LHC experiments soon. At the LHC, running at 7 TeV, with an integrated luminosity of only ∼ 200 pb −1 , it should be possible to record about 30000 top quark pair events before selection [4]. With the planned increases in center of mass energy and luminosity, the LHC will turn into a veritable top-quark factory, producing millions of top-quark events per year [5].
In order to take full advantage of the improved experimental measurements, it will be necessary to provide theoretical predictions for the measured observables which are as accurate as possible. Most of the observables related to the top-quark pair production have been calculated up to NLO [6]. In several cases, the next-to-leading (NLO) corrections have been supplemented by the resummation of large logarithmic corrections at leading (LL, [7]), next-to-leading (NLL, [8]) and next-to-next-to-leading logarithmic (NNLL, [9][10][11][12]) accuracy. However, to match the precision of the forthcoming experimental data, full next-to-next-to-leading order (NNLO) calculations are required for at least some of the observables, such as the top-quark pair production total cross section [13].
To obtain NNLO predictions for the top quark pair production in QCD, the following matrix elements need to be computed: i) two-loop matrix elements in both the quarkantiquark annihilation channel and the gluon-fusion channel; ii) one-loop matrix elements with one extra parton in the final state; and iii) tree-level matrix elements with two extra partons in the final state. The diagrams belonging to the sets ii) and iii) form part of the NLO corrections to top-pair-plus-jet production and have been calculated by several groups [14][15][16]. Aside from the calculation of the required Feynman diagrams, a full NNLO calculation of the top-quark pair production cross section is particularly challenging because it entails the development of a NNLO subtraction method, which could be used in processes with massive partons [17][18][19][20][21].
As far as matrix elements are concerned, the last missing ingredient to the NNLO corrections to top quark pair production are the two-loop matrix elements for qq → tt and gg → tt. Both types of matrix elements were calculated in the s ≫ m 2 limit (where s is the partonic center of mass energy and m is the mass of the top quark) [22,23]. A fully numerical calculation of the two-loop corrections in the quark-antiquark annihilation channel is also available [24]. This paper is the third of a series of works aiming towards the analytic calculation of the two-loop matrix elements for the top pair production process. For what concerns the quark-antiquark annihilation channel, the two-loop diagrams involving a closed light or heavy-quark loop were evaluated in [25], while the two-loop diagrams contributing to the leading color coefficient were evaluated in [26]. In both cases, the results obtained retain the full dependence on the top-quark mass and on the kinematic invariants; they agree with the numerical results of [24]. Having analytical results available has several advantages over a purely numerical representation. Besides their considerably shorter evaluation time, the analytical results also allow for an expansion in different kinematical limits (threshold, high energy). In the present paper, an analytical expression for the two-loop diagrams contributing to the leading color coefficient in the gluon-fusion channel is derived. We carry out the calculation by employing the technique based on the Laporta algorithm [27] and the differential equation method [28], already used in [25,26]. The calculation of the leading color coefficient in the gluon fusion does not require the calculation of any new master integrals beyond the ones obtained in the two previous works, such that we do not discuss the calculational method in full detail. The interested reader can find in [25,26] a detailed description of the techniques employed.
The paper is organized as follows: in Section 2, we introduce our notation and conventions; in Section 3, we summarize the most relevant features of our calculational method. Section 4 describes the UV renormalization of the bare amplitude. The resulting two-loop amplitude contributions are described in Section 5, where we also provide numerical values in some benchmark points, and discuss the expansion in the threshold limit. Finally, Section 6 contains our conclusions.

Notation and Conventions
We consider the scattering process in Euclidean kinematics, where p 2 i = 0 for i = 1, 2 and p 2 j = −m 2 for j = 3, 4. The Mandelstam variables are defined as follows (2.2) Conservation of momentum implies that s + t + u = 2m 2 . The squared matrix element (summed over spin and color), calculated in d = 4 − 2ε dimensions, can be expanded in powers of the strong coupling constant α S as follows: The tree-level amplitude involves the three diagrams shown in Fig. 1 and their contribution to Eq. (2.3) is given by 4) where N c is the number of colors, 3) arises from the interference of one-loop diagrams with the tree-level amplitude [6]. The O(α 2 S ) term A 2 consists of two parts, the interference of two-loop diagrams with the Born amplitude and the interference of one-loop diagrams among themselves: was derived in [29][30][31][32]. A (2×0) 2 , originating from the two-loop diagrams, can be decomposed according to color and flavor structures as follows: 6) where N l and N h are the number of light-and heavy-quark flavors, respectively. The coefficients A, B, . . . , I lh in Eq. (2.6) are functions of s, t, and m, as well as of the dimensional regulator ε. These quantities were calculated in [23] in the approximation s, |t|, |u| ≫ m 2 . For a fully differential description of top quark pair production at NNLO, the complete mass dependence of A (2×0) 2 is required; to date no such a result is available. Recently, a formula which allows to calculate the IR poles of a generic two-loop QCD amplitude was derived [33,34] and employed to obtain analytic expressions for all of the IR poles for the top-quark pair production [35].
In this work, we provide an exact analytic expression for the leading colour coefficient A in Eq. (2.6), which receives contributions from planar Feynman diagrams only.

Calculation
The package QGRAF [36] generates 789 two-loop Feynman diagrams contributing to the process gg → tt (considering one massless and one massive flavor). Evaluating their color structures, we find that 300 of them contribute to the leading color coefficient in Eq. (2.6). Since we use the covariant sum over the polarizations of the incoming gluons, we have to consider additional 116 (out of 218) diagrams for the process initiated by ghosts. We interfere the two-loop diagrams with the tree-level amplitude, and calculate the color and Dirac traces before carrying out the integration over the loop momenta. The resulting scalar loop integrals are reduced to a set of Master Integrals (MIs) employing the technique based upon the Laporta algorithm [27]. Then, the MIs can be evaluated analytically by means of the differential equation method [28].
Starting with the QGRAF output, the calculation of the interferences in terms of the MIs was carried out with a parallelized C++ package: Reduze 2 [37]. In particular, this code provides a fully distributed variant of the Laporta algorithm for the reduction 1 of the loop integrals. The package employs GiNaC [41] for the algebraic manipulations.
The MIs needed for the calculation presented in this paper were already known in the literature [25,26,[42][43][44][45][46][47][48]. In particular, all of the two-loop four-point MIs encountered in the calculation of the leading color coefficient in the gluon-fusion channel coincide with the ones needed for the corresponding calculation in the quark-antiquark annihilation channel [26], or can be obtained by the latter by replacing the Mandelstam variable t with u. All of the MIs were calculated in the non-physical region s < 0, where they are real. The result in the physical region is then recovered by analytic continuation.
The transcendental functions appearing in the results are one-and two-dimensional harmonic polylogarithms (HPLs) [46,[49][50][51] of maximum weight four. All of the HPLs appearing in the analytic expression of the coefficient A can be evaluated numerically with arbitrary precision by employing the methods and codes described in [50]. A detailed list of the functional basis employed in this calculation can be found in Appendix A of [26]. Appendix B in the same paper describes the expansion of this class of HPLs in the threshold limit β → 0, which requires some care. As expected, the sum of the bare two-loop corrections (as well as the UV renormalized corrections) is symmetric with respect to the exchange t ↔ u.

Renormalization
The renormalized QCD amplitude is obtained from the bare one by expanding the following 1 Other reduction codes released publicly are the Maple package A.I.R. [38], the Mathematica package FIRE [39] and the C++ package Reduze [40].

expression :
where Z WF,n (n = g, t) are the external particle wave-function renormalization factors, α S is the renormalized coupling constant, Z α S is the coupling constant renormalization factor, m is the renormalized heavy-quark mass, and Z m is the mass renormalization factor. In the rest of the section we suppress the subscript "S" in α S . We introduce the following auxiliary quantities: By expanding the amplitude and the wave function renormalization factor in a 0 we find: The relation between a 0 and a is given by: By employing Eqs. (4.4,4.5) in Eq. (4.2) we find ren + O(a 4 ) , A (1) In the equations above, A i , i = 0, 1, 2, represents the bare i-loop amplitude stripped of the factor a, A where C(ε) = (4π) ε Γ(1 + ε), β 0 = (11/6)C A − (1/3)(N l + N h ) and where γ is the Euler-Mascheroni constant γ ≈ 0.577216. Concerning the two-loop renormalization factors, we provide only the terms which contribute to the renormalization of the leading color coefficient. They are (see for instance [52]): (4.11) δZ (2) WF,g = 0 , (4.12) . (4.14)

Results
The main result of the present paper is an analytic expression for the coefficient A in Eq. (2.6) which retains the complete dependence on the variables t, u, on the renormalization scale µ, and on the top-quark mass m. The analytic result is too long to be explicitly presented in written form. To make it available to the reader we include in the arXiv submission of this work a text file with the complete result. The coefficients in Eq. (2.6) still contain infrared poles. This makes our result dependent on the choice of a global, ε-dependent normalization factor. We choose to present the coefficient A factoring out an overall term Another possible normalization, widely used in the literature, is the one in which a term is retained as an overall factor (see for instance [35] and [24]). Our result, therefore, can be compared to the results in [24,35] by multiplying it by a coefficient While the formulae in the text and the supplied files are given in the normalization defined by Eq. (5.1), for ease of comparison with other results in the literature we present in Table 1 and in Figures 3-4 a few benchmark points and plots in the normalization of Eq. (5.2). We provide a code which numerically evaluates the analytic expression of the leading color coefficient for arbitrary values of the mass scales involved in the calculation. The code is written in C++ and uses the package for the evaluation of multiple polylogarithms within GiNaC [50].
Our result was cross checked by comparing it with the partial results already available in the literature. By expanding the analytic expression of the coefficient A in the limit s, |t|, |u| ≫ m 2 and neglecting terms suppressed by positive powers of m 2 /s for t/s = const , we obtained the result published in [23]. For the IR poles of the coefficient A we find numerical agreement with the results of [35].
In Fig. 3 the finite part of A is plotted as a surface depending on the variables η and φ, defined as Finally, it is possible to expand the expression for A for values of the center of mass energy close to the production threshold. We define where θ is the scattering angle in the partonic center of mass frame. Keeping ξ = const , we expand our results in powers of the heavy quark velocity β up to terms of order β 5 included. The coefficients of this expansion contain transcendental constants which originate from one-and two-dimensional HPLs evaluated at x = 1. Since we did not find a satisfactory analytical representation for all of these constants, in the formulae below we provide the coefficients of the expansion in a numerical form. We find: with   We observe that the dependence on β and on ξ in the formulas above is polynomial. Moreover, β enters only via powers of β 2 and ξ only via powers of ξ(1 − ξ). The latter dependence explicitely reflects the symmetry of A under exchange of forward and backward directions, ξ → 1 − ξ. All of the logarithmic terms ln β, ln ξ, ln (1 − ξ), ln (1 − 2 ξ), . . ., which are indeed present in the expansion of individual HPLs, cancel out in the final expressions. Thus, the coefficient A is finite at threshold. The expansion presented here could be used in the future for the calculation of logarithmically enhanced leading colour terms near the tt production threshold.

9)
In Fig. 4, we show the comparison between the exact expression of the coefficient A and the expansions in the two regimes: threshold expansion on the left hand side and small-mass expansion on the right hand side. For the small-mass expansion we consider the limit m 2 /s → 0 for φ = −(t − m 2 )/s = const and neglect terms suppressed by powers of m 2 /s larger than a given order. This prescription preserves the symmetry of the result under exchange of forward and backward directions, φ → 1 − φ, for a given order of the expansion.
With the arXiv submission of this paper, we provide algebraic expressions and numerical implementations of the threshold expansion up to order β 5 included and of the small mass expansion up to order (m 2 /s) 2 included.

Conclusions and Outlook
In this paper, we presented the analytic calculation of the two-loop leading color corrections to the heavy-quark production matrix element in the gluon fusion channel. The diagrams required to calculate this coefficient are all planar. The result presented here retains the exact heavy-quark mass dependence; no assumptions on the hierarchy between the mass scales involved in the problem was made. The formula for the coefficient A in Eq. (2.6), which was obtained in this work, was validated against analytic results valid in the small-mass region [23]. In addition, we found numerical agreement with the exact analytic expression for the IR poles which was obtained in [35].
Our result represents a gauge invariant sub-set of the full two-loop corrections to the partonic process gg → tt. In order to complete the analytic calculation of the twoloop corrections, it is necessary to calculate all of the fermionic diagrams, and the nonplanar gluonic diagrams. A large part of the non-leading color coefficients in Eq. (2.6) can be calculated within the same calculational framework employed here. However, it is known that some of the two-loop diagrams appearing in the gluon fusion channel cannot be expressed in terms the HPLs functional basis. For example, some of the diagrams with a closed heavy-quark loop involve a "sunrise"-type subtopology with three equal massive propagators and an external momentum which is not on the mass shell of the internal propagators. Such a three-propagator graph can be written only in terms of elliptic integrals [53]. The use of numerical [24] or seminumerical [54,55] methods could thus be unavoidable in the evaluation of these diagrams.
In order to obtain NNLO predictions for the total tt production cross section and for differential distributions, it is necessary to combine the two-loop virtual corrections with the already available one-loop corrections to the tt+(1 parton) process and with the tree-level matrix elements for the process tt+(2 partons) [14][15][16]. These diagrams with additional partons in the final state contribute to infrared-divergent configurations where up to two partons can become unresolved. Their implementation requires the application of a NNLO subtraction method. The methods presently available [17][18][19] do not provide yet the full counterterms for a 2 → 2 hadronic process involving massive partons. Two subtraction methods for the NNLO calculation of the top-quark pair production cross section were outlined recently [20,21].
In the light of these advancements, the calculation of the full NNLO corrections to top quark pair production in hadronic collisions is gradually becoming feasible.