Heavy-quark production in gluon fusion at two loops in QCD

We present the two-loop virtual QCD corrections to the production of heavy quarks in gluon fusion. The results are exact in the limit when all kinematical invariants are large compared to the mass of the heavy quark up to terms suppressed by powers of the heavy-quark mass. Our derivation uses a simple relation between massless and massive QCD scattering amplitudes as well as a direct calculation of the massive amplitude at two loops. The results presented here together with those obtained previously for quark-quark scattering form important parts of the next-to-next-to-leading order QCD corrections to heavy-quark production in hadron-hadron collisions.


Introduction
The production of heavy quarks at hadron colliders is an important process for a range of reasons. Experimentally, events with a top-quark, being the heaviest quark known thus far, lead to very characteristic signatures in a hadronic scattering process. This allows for event reconstruction in a variety of channels, which e.g. at LHC is supplemented by an anticipated very high statistics for the production cross section [1,2]. In this way precision measurements of the top-mass, of (differential) tt-distributions and also tests of the production and the subsequent decay mechanism (including anomalous couplings and top-spin correlations) can be conducted at Tevatron or LHC. The precise knowledge of the top-quark parameters has direct impact on precision tests of the Standard Model, in particular on the Higgs sector due to the large top-Yukawa coupling. Moreover, events with top-quarks will also make up for a large part of the background in searches for the Higgs boson or new physics.
In the case of bottom-quark production important measurements are inclusive B-meson and bjet production at moderate to large transverse momentum. At colliders b-flavored jets are produced for instance in decays of top-quarks, the Higgs boson and numerous particles proposed in extensions of the Standard Model. Differential b-jet distributions are also important for measurements of parton distributions. Predictions for heavy-quark hadro-production have theoretical uncertainties, which can even be larger than e.g. the corresponding predictions for the production of light jets. This is due to radiative corrections in Quantum Chromodynamics (QCD), where typically large logarithms at higher orders appear. At present, the next-to-leading order (NLO) QCD corrections for heavy-quark hadro-production are known [3][4][5][6][7][8][9] and, once they become available, the complete next-to-next-to-leading order (NNLO) QCD corrections will reduce the uncertainty of theory predictions by improving the stability with respect to scale variations. For inclusive B-meson production [10,11], for instance, the formalism of perturbative fragmentation functions of a heavy quark [12] has already been extended through NNLO by providing the initial conditions [13,14] together with the time-like (non-singlet) evolution [15]. Alternative ways devoted to control the theoretical uncertainty in the case of b-quarks, for example through the definition of dedicated jet algorithms have been proposed recently [16].
In the hadronic production of heavy quarks large perturbative QCD corrections arise in different kinematical regions. If the heavy quarks are produced with (partonic) center-of-mass energies s close to threshold, then Sudakov-type logarithms, powers of log(β), appear to all orders in perturbation theory, where β = 1 − 4m 2 /s is the velocity of the heavy quark. The necessary resummation has been carried out to next-to-leading logarithmic (NLL) accuracy [17][18][19] and has successfully improved the phenomenology of top-quark production at Tevatron. In the high-energy regime on the other hand, when the heavy quarks are fast, logarithms in the quark mass m appear. These log(m)-terms are of collinear origin and dominate cross sections when m becomes negligible in comparison to other kinematical invariants, such as e.g. for the inclusive b-jet spectrum at large transverse momentum. All-order predictions can be achieved by means of an explicit resummation. The necessary technology for resumming both the incoming and outgoing collinear logarithms to NLL accuracy is available, see e.g. Refs. [20,21] and references therein.
In this article, we want to further improve the precision of fixed-order perturbative QCD corrections. To that end, we present results for the virtual QCD corrections at two loops for the pair-production of heavy quarks in gluon fusion. Together with the corresponding results for quark-quark scattering obtained previously by us [22] these are essential parts of the complete NNLO QCD corrections. To be precise, we calculate the interference of the two-loop with the Born amplitude. We work in the limit of fixed scattering angle and high energy, where all kinematic invariants are large compared to the heavy-quark mass m. Thus, our result contains all logarithms log(m) as well as all constant contributions (i.e. the mass-independent terms) and we consistently neglect power corrections in m.
In our calculation we employ two different methods. On the one hand, we apply a generalization of the infrared factorization formula for massless QCD amplitudes [23,24] to the case of massive partons [25]. In a nut-shell this results in an extremely simple universal multiplicative relation between a massive QCD amplitude in the small-mass limit and its massless version [25]. In this way, we can largely use for our derivation the results of the NNLO QCD corrections to massless quark-gluon scattering (i.e. gg → qq). These have been computed for the squared matrix element, i.e. the interference of the two-loop virtual corrections with the Born amplitude in Ref. [26], while results for the individual (independent) helicity configurations of the two-loop amplitude |M (2) itself have been given in Refs. [27,28]. On the other hand, we perform a direct calculation of the relevant Feynman diagrams in the massive case followed by a subsequent expansion in the small-mass limit by means of Mellin-Barnes representations.
The outline of the article is as follows. In Section 2 we give some basic definitions and present a short summary of our methods in Section 3. There we briefly explain the essence of the QCD factorization approach to calculate massive amplitudes and give the relevant formulae. We also highlight the key steps of the direct evaluation of Feynman diagrams in the small-mass limit and, in particular, comment on the non-planar topologies. Section 4 contains our results and we conclude in Section 5. Appendix A gives the explicit result for a massive non-planar scalar integral in the small mass limit.

Setting the stage
The pair-production of heavy quarks in the gluon fusion process corresponds to the scattering where p i denote the gluon and quark momenta and m the mass of the heavy quark. Energymomentum conservation implies p Following the notation of Ref. [26] we consider the scattering amplitude M for the process (1) at fixed values of the external parton momenta p i , thus p 2 1 = p 2 2 = 0 and p 2 3 = p 2 4 = m 2 . The amplitude M may be written as a series expansion in the strong coupling α s , and we define the expansion coefficients in powers of α s (µ 2 )/(2π) with µ being the renormalization scale. We work in conventional dimensional regularization, d = 4 − 2ε, in the MS-scheme for the coupling constant renormalization. The heavy mass m on the other hand is always taken to be the pole mass.
We explicitly relate the bare (unrenormalized) coupling α b s to the renormalized coupling α s by where we put the factor S ε = (4π) ε exp(−εγ E ) = 1 for simplicity and β is the QCD β-function [29,30] The color factors are in a non-Abelian SU(N)-gauge theory C A = N, C F = (N 2 − 1)/2N and T F = 1/2. Throughout this letter, N denotes the number of colors and n f the total number of flavors, which is the sum of n l light and n h heavy quarks.
In the following, we will confine ourselves to the discussion of the squared amplitude for the process (1), although it should be clear that our approach and the results of the present article can be easily extended to the (color ordered) partial amplitudes for the individual helicity combinations of the massive two-loop amplitude |M (2) itself. There one would rely in particular on Refs. [27,28].

(6)
A is a function of the Mandelstam variables s, t and u given by and has a perturbative expansion similar to Eq. (3), In terms of the amplitudes the expansion coefficients in Eq. (8) may be expressed as where we have discarded powers in the heavy-quark mass m in A 4 .
The expressions for A 6 have been presented e.g. in Refs. [8,9] and the loop-by-loop contribution M (1) |M (1) at NNLO in A 8 has been published in Ref. [31]. Both results for A 6 and M (1) |M (1) have been obtained in dimensional regularization and with the complete dependence on the heavy-quark mass. Here we provide for the first time the real part of M (0) |M (2) up to powers O (m) in the heavy-quark mass m.

The massive amplitude from QCD factorization
Let us briefly recall the key findings of Ref. [25] on how to calculate loop amplitudes with massive partons from purely massless amplitudes. Heuristically, the QCD factorization approach rests on the fact that a massive amplitude M (m) for any given physical process shares essential properties in the small-mass limit with the corresponding massless amplitude M (m=0) . The latter one, M (m=0) , generally displays two types of singularities, soft and collinear, related to the emission of gluons with vanishing energy and to collinear parton radiation off massless hard partons, respectively. These appear explicitly as poles in ε in dimensional regularization after the usual ultraviolet renormalization is performed. In the former case, the soft singularities remain in M (m) as single poles in ε while some of the collinear singularities are now screened by the mass m of the heavy fields, which gives rise to a logarithmic dependence on m, see e.g. Ref. [32].
This structure of singularities for massless amplitudes has been clarified to all orders in perturbation theory [24] as all 1/ε terms can be exponentiated, see also Ref. [33]. Similarly, all poles in ε and log(m) terms for amplitudes with massive partons also obey an all-order exponentiation with mostly the same anomalous dimensions as in the massless case [25]. Thus, in the small-mass limit the differences between a massless and a massive amplitude can be thought of as due to the difference in the infrared regularization schemes. QCD factorization provides a remarkably simple direct relation between M (m) and M (m=0) The function Z (m|0) is process independent and depends only on the type of external parton, i.e. quarks and gluons in the case at hand. For external massive quarks Q it is defined as the ratio of the on-shell heavy-quark form factor and the massless on-shell one, both being known [34][35][36] to sufficient orders in α s and powers of ε. An explicit expression for up to two loops is given in Ref. [25] (note the different normalization α s /(4π) used there). The leading n f terms ∼ (n f α s ) n for our process Eq. (1), gg → QQ, can also be predicted based on the above arguments. Keeping only terms quadratic in n h and/or n f = n h + n l one has up to two loops: where Z (2) [g] = Z with Z [g] ∼ n h given in Ref. [25] (again note the different normalization α s /(4π) used there). Z (1) [g] is also equal to the O (α s ) term in the gluon wave function renormalization constant Z 3 in Eq. (22).
to Z 3 is discussed after Eq. (24) below. To derive Eq. (15) we apply the definition for Z (m|0) given in Ref. [25], i.e. evaluate the ratio of the gluon form factor with heavyloop insertions and the pure massless gluon form factor [36,37]. The additional renormalization constant that enters the effective Hgg vertex (see e.g. Ref. [37] for details) cancels in the ratio and does not contribute to Z (m|0) [g] . Exploiting the predictive power of the relation Eq. (12) and applying it to the process Eq. (1) we get (16) which assumes the hierarchy of scales m 2 ≪ s,t, u , i.e. we neglect terms O (m). Eq. (16) predicts the complete real part of the squared amplitude M (0) |M (2) (m) except (as indicated) for those terms, which are linear in n h (the number of heavy quarks) and, at the same time not proportional to n l (the number of light quarks). These two-loop contributions have been excluded explicitly also from the definition [25] of Z (m|0) , as one needs additional process dependent terms for their description.
The two-loop massless amplitudes Re M (0) |M (2) (m=0) are computed in Ref. [26]. We have checked that the finite remainders of the squared two-loop amplitudes obtained after the infrared subtraction procedure discussed in that reference agree with the corresponding terms constructed from the two-loop helicity amplitudes calculated in Ref. [27]. We have also found similar agreement between the finite remainders of the qq → q ′q′ amplitudes we used in Ref. [22] that we extracted from Refs. [38,39].

Direct computation of the massive amplitude
The direct computation of the massive amplitude proceeds according to the same scheme as in our previous publication [22], which itself was an evolution of the methodology developed in Refs. [40][41][42]. In short, the complete amplitude is reduced to an expression containing only a small number of integrals with the help of the Laporta algorithm [43]. In a next step, Mellin-Barnes (MB) representations [44,45] of all these integrals are constructed, and analytically continued in the dimension of space-time with the help of the MB package [46] revealing the full singularity structure. A subsequent asymptotic expansion in the mass parameter is done by closing contours and resumming the integrals, either with the help of XSummer [47], or the PSLQ algorithm [48].
We shall now concentrate on the differences with respect to our previous calculation. First of all, the number of master integrals is substantially larger, reaching 422, which adds a lot to the complexity of the computation. This is partly due to the fact, that the symmetry with respect to the exchange of the gluons generates the same topologies in the t-and u-channels, but more importantly because of completely new topologies, which come together with the more extended set of gluon interactions. Whereas in Ref. [22], it was possible to avoid the computation of the high-energy asymptotics of non-planar graphs, and still have a test of the factorization approach, we were not able to avoid them here. In fact, this additional complication is due to the single heavyquark loop diagrams of Fig. 1, which are explicitely removed from the factorization approach. The complete set of non-planar master integrals belonging to this class is depicted in Fig. 3.  The first problem that one has to face when dealing with non-planar integrals is the construction of MB representations. In the planar case, the iterative loop-by-loop integration has proved to be the most fruitful. On the other hand, the first non-planar double-box diagram ever computed [45], with massless propagators and on-shell external legs, had its four-dimensional MB representation, I m=0, dim=4 NP , derived directly from the two-loop Feynman parameter representation. It seems, however, that when masses are involved, the loop-by-loop representations are more compact, as seen for example in Ref. [49,50]. However, any asymptotic expansion contains a so-called hard part, which is obtained by setting all the small parameters to zero, and which in our case would correspond to the massless graph. Following this line of thought, one can derive a representation for the massless on-shell graph by taking suitable residues in the result presented in Ref. [49]. Unfortunately, one arrives at a six-fold representation, I m=0, dim=6 NP , in clear disadvantage with respect to I m=0, dim=4 NP . Interestingly, there is an even more severe problem inherent in the loop-by-loop approach. The leading pole derived in Ref. [45] reads (up to normalization factors irrelevant for this discussion) whereas the six-fold representation gives (with the same normalization) This clear discrepancy is only explained when we look at the subleading pole from the six-fold representation, which contains a logarithmic singularity Obviously, the extension of the integral into the Euclidean domain, performed with the help of the u parameter, regularized part of the infrared singularity. The only way to obtain the correct result would be to first take the limit u → −s − t, and only then ε → 0. This, however, is a highly non-trivial task, as is well known from studies aiming at the derivation of exact expressions of Feynman integrals in d-dimensions.
In view of all the above arguments, we derived our MB representations directly from the twoloop Feynman parametric representation. In particular, for the scalar integral of Fig. 2, we have the following six-fold representation where the loop integration is done with the measure e εγ E R d d k/(iπ d/2 ) per loop. We defer the presentation of the full result of the expansion of this integral to the Appendix A. Here we only note, that the leading term of the expansion has a square root singularity, which is a feature of non-planar graphs, that does not occur in any of the planar integrals considered in this calculation. Clearly, the disappearance of this square root singularity is a simple test of the correctness of the calculation.
The second problem that requires care is connected with the choice of the master integrals. In Fig. 3, we have not only shown the topologies, but also the multiplicities of the masters. The basic seven-liner needs as much as five different integrals. It is clear that we want to avoid coefficients containing poles in ε or m 2 , since the leading behavior of the integrals is difficult enough to determine, and such poles would be synonymous of higher orders in the respective expansions. After inspection it turns out that one can take two integrals with second powers of the denominators into the set, but tensors rank one and two are unavoidable. Although it is possible to generate representations for arbitrary tensors, one ends up with a huge number of four-fold integrals after expansion. Instead, one can introduce a new ficticious propagator that will have a negative power. For this we choose the square of the momentum that runs through the crossed box subloop, as shown in Fig. 2. The final set of seven-line non-planar master integrals is shown in Fig. 4. Note that the dotted masters, i.e. those having higher powers of chosen propagators are particularly easy to calculate, because one only needs the singularities in 1/m 2 . All of the representations can be derived from the following one where a = ∑ 8 i=1 a i , a S = ∑ i∈S a i with S a subset of 1, ..., 8, and a i with i = 1, ..., 7 are the powers of the denominators according to the labeling given in Fig. 2, whereas a 8 is the power of the additional denominator 1/k 2 . Even though Eq. (21) has two more integrations than Eq. (20), the presence of the factor 1/Γ(a 8 ) makes it necessary to perform first an analytic continuation in a 8 to a negative integer value, which effectively reduces the number of integration variables back to six. In fact, the above representation can be used to compute any of the master integrals from the full set of non-planars of Fig. 3.
To complete our exposition of the direct calculation of the amplitude, we have to note that the renormalization of the bare amplitude requires the on-shell wave function renormalization constant of the gluon, which is non-vanishing due to the presence of heavy-quark loops. We give it here as an expression exact in d-dimensions, where a b s is defined by the bare coupling constant α b s and with S ε = (4π) ε exp(−εγ E ), and In terms of the renormalized coupling α s and expanding in powers of ε through sufficient terms, the result for Z 3 reads The corresponding expression for the wave function renormalization constant Z 2 of a light quark state has been given in Ref. [22].
Our result for Z 3 in Eqs. (22) and (24) coincides to first order in α s with Z [g] from Eq. (14). At second order in α s all known terms (i.e. those quadratic in the number of flavors) in the two constants are also a complete match (also in d-dimensions). In our direct evaluation of Z 3 in Eq. (22) we observe gauge independence through two loops (within the class of covariant gauges employed in the calculation), which is consistent with the arguments given in Ref. [25] in favor of the identification of the constant Z (m|0) [g] with Z 3 evaluated in a physical gauge.

Results
We are now ready to present our result for gg → QQ scattering for the interference of the two-loop and Born amplitude, which we have ordered according to the power of the number of colors N and the numbers of n l light and n h = 1 heavy quarks with n f = n l + n h total flavors.
The coefficients A, E l , H l , H lh , H h , I l , I lh , and I h have been computed with both of our methods. We have found agreement between the direct computation of the relevant Feynman diagrams in the small-mass expansion and the QCD factorization approach as given by the universal multiplicative relation (12). All other terms linear in n h , that includes E h , F h and G h have been obtained by means of a direct calculations of the massive loop integrals as detailed above. The remaining coefficients B, C, D, F l and G l have been derived by application of the factorization formula.
We choose x = −t/s as the only dimensionless kinematic variable in the problem and we keep the dependence on the renormalization scale, µ, explicit. We also introduce the following compact notation The different components now read consistent NNLO treatment. At one-loop the complete structure of the singularities as well as all large logarithms in the heavy-quark mass m can therefore be entirely treated with the methods of Ref. [32].
Our derivation relies on the combination of two completely different methods and we have ensured substantial overlap between them. In this way we have had mutual and highly non-trivial cross checks on our direct calculation of massive Feynman diagrams and on the QCD factorization approach [25]. In particular through the latter method we have had the possibility to relate to various massless [26-28, 35, 36] and massive results [34] available in the literature and we have found consistency.
The results for gluon fusion in the present article (and Ref. [22] for quark-quark scattering) are of direct relevance whenever power corrections in the heavy-quark mass are negligible. This is certainly the case for hadro-production of bottom pairs over a large kinematical range at colliders and, to a lesser extent perhaps, for tt-production at LHC energies. Here, possible improvements would come from the systematic computation of power corrections in the heavy-quark mass, i.e. terms proportional to (m 2 /s) k with k ≥ 1 to improve the convergence of the small-mass expansion. This could be achieved, for instance, by extending the methods of Sec. 3 for the direct calculation of massive Feynman diagrams to higher powers in m.
Finally, it is clear, that our result for M (0) |M (2) still has to be combined with the tree-level 2 → 4, the one-loop 2 → 3 as well as the square of the one-loop 2 → 2 processes M (1) |M (1) in order to yield physical cross sections. Some of the matrix elements including the full mass dependence can be easily generated, others have become available in the literature only rather recently, see e.g. Refs. [31,51]. The combination of all these contributions enables the analytic cancellation of the remaining infrared divergences as well as the isolation of the initial state singularities. The latter will have to be absorbed into parton distribution functions of, say, the proton in order to match with a precise parton evolution at NNLO [52,53]. All these remaining steps are necessary prerequisites e.g. to the construction of numerical programs which then provide NNLO QCD estimates of observable scattering cross sections for heavy-quark hadro-production.
A MATHEMATICA file with our results can be obtained by downloading the source from the preprint server http://arXiv.org. The results are also available from the authors upon request.

A Asymptotics of the massive non-planar scalar integral
Here, we give the leading high energy behavior of the non-planar scalar integral with a massive loop corresponding to the MB representation Eq. (20)