Two-loop scattering amplitude for heavy-quark pair production through light-quark annihilation in QCD

We present the first full analytic evaluation of the scattering amplitude for the process $q {\bar q} \to Q {\bar Q}$ up-to two loops in Quantum Chromodynamics, for a massless $(q)$ and a massive $(Q)$ quark flavour. The interference terms of the one- and two-loop amplitudes with the Born amplitude, decomposed in terms of gauge invariant form factors depending on the colour and flavour structure, are analytically calculated by keeping complete dependence on the squared center-of-mass energy, the squared momentum transfer, and the heavy-quark mass. The results are expressed as Laurent series around four space-time dimensions, with coefficients given in terms of generalised polylogarithms and transcendental constants up-to weight four. Our results validate the known, purely numerical calculations of the squared amplitude, and extend the analytic knowledge, previously limited to a subset of form factors, to their whole set, coming from both planar and non-planar diagrams, up-to the second order corrections in the strong coupling constant.


Introduction
The production of top quark (t) and its anti-particle (t) has been occupying a central role within the precision physics programme at hadron colliders, like the Tevatron and the Large Hadron Collider (LHC), over the last three decades. Being the heaviest known elementary particle, the t quark has offered a portal to the discovery of the Higgs boson, and it is considered pivotal for understanding the electroweak symmetry breaking mechanism. Studies of top-quark production (and decay) at the current LHC physics programme enters the high-precision tests of the parameters of the Standard Model (SM), such as couplings and masses, as well as the analyses of backgrounds, for discriminating deviations that could indicate the path to move beyond it. Within SM, the production of tt pairs in hadronic collisions is the main source of top quarks, therefore, it is considered among the cornerstone processes at the current and future hadron colliders. Because of its role for the precision physics programme at hadron colliders, the tt-pair production has triggered a significant progress in the developments of theoretical methods for determination of the (differential) cross-sections, hence it has been stimulating the constant effort of providing calculations in Perturbative Quantum Chromodynamics (QCD), of increasing order in the strong-coupling series expansion.
The calculation of the NNLO QCD corrections to pp → tt requires four types of terms: the double-real corrections, coming from the tree-level squared amplitude for a process with two additional partons in the final state; the real-virtual corrections, due to the interference of the tree-level and of the one-loop amplitude for a process with one additional gluon in the final state; the squared one-loop corrections; the double-virtual corrections, due to the interference of the two-loop amplitude with the tree-level one.
The scattering amplitude for the real-virtual contributions were evaluated in [30,31], and more recently in [32]. The purely virtual contributions depend on the square of one-loop amplitude and the genuine two-loop amplitude. The former has been computed analytically in [33][34][35], while the latter has been determined completely numerically in [36][37][38]. The analytic evaluation of the two-loop amplitude is known partially [39][40][41][42][43][44]. The main difficulty, in this case, is due to the analytic evaluation of the independent integrals appearing in the decomposition of the two-loop amplitudes, known as master integrals (MIs).
At parton-level, the tt-production proceeds via the annihilation of a light-quark (q) and an anti-quark (q), qq → tt, and the more luminous gluon-fusion channel, gg → tt. As regarding the gluon-fusion channel, the analytic evaluation of the interference of the two-loop amplitude with the tree-level amplitude is only partially complete, and they are expressed in terms of generalised polylogarithms (GPLs) and elliptic integrals [41][42][43][45][46][47]. Very recently, the two-loop helicity amplitudes for the tt-production in the gluon-fusion channel within the leading colour approximation, including the contribution of closed loops of quarks, has been computed in [44]. As regarding the light-quark pair annihilation channel, the interference of the two-loop amplitude with the corresponding tree-level amplitude can be decomposed in terms of ten form factors, according to the colour and flavour structure. Eight of them are known analytically, and expressed in terms of GPLs [39,40,48].
In this work, we present the complete analytic evaluation of the two-loop scattering amplitude for the scattering process qq → QQ, with a massless (q) and a massive (Q) quark flavour, in QCD, including leading and sub-leading colour contributions. We calculate the whole set of ten form factors analytically, including the two form factors previously unavailable, which take contribution from both planar and non-planar graphs. The latter do not contribute to the eight form factors already known, and their evaluation constitute part of the novel insights of the current work.
The loop integrals appearing in the un-renormalised interference terms of the oneand two-loop bare amplitudes with the leading-order one are regulated within the Conventional Dimensional Regularisation (CDR), where d is the number of continuous space-time dimensions.
The calculation is automated within the Aida [49] framework, implementing the adaptive integrand decomposition algorithm [50,51] and interfaced: to FeynArts [52], FeynCalc [53], for the automatic diagram generation and algebraic manipulations of the integrands; to Reduze [54], and Kira [55], for the generation of the relations required for the decomposition in terms of MIs; to SecDec [56], for the numerical evaluation of MIs, if needed; to PolylogTools [57], Ginac [58], and HandyG [59], for the numerical evaluation of the analytic expressions. The cancellation of the ultraviolet (UV) divergences of the bare interference terms at one and two loops is carried out by renormalising the quark fields and masses in the on-shell scheme, and the strong coupling in the MS-scheme, along the lines of [36,39]. By using the analytic expressions of the MIs [39,[60][61][62][63][64][65], the renormalised interference terms are finally expressed as Laurent series around d = 4 dimensions, by keeping the complete dependence on the Mandelstam invariants s and t, and on the heavyquark mass M . The one-and two-loop contributions are computed, respectively, up-to the first-order term, and up-to the finite term, in the four dimensional series expansion, whose coefficients are expressed in terms of GPLs and transcendental constants of up-to weight four. The analytic results are obtained in a non-physical region, where the variables s and t are negative, and are numerically continued to the physical region, above the heavy-quark pair-production threshold, s ≥ 4M 2 .
The study of the virtual NNLO QCD corrections for the process qq → QQ, hereby presented, extends to the non-Abelian case the study of the four-fermion scattering amplitude with one massive fermion pair, in gauge theories, recently completed for the process e + e − → µ + µ − in Quantum Electrodynamics (QED) [79].
In the following pages, we describe the strategy we adopted to solve the problem of the analytic evaluation of the double-virtual NNLO corrections to one, out of two, partonic reactions contributing the hadroproduction of heavy-quark pair. Thus, providing what we consider an important validation and extension of the purely numerically known results, which have been employed to obtain state-of-the-art perturbative predictions within topquark physics studies at hadron colliders (see [80,81] and reference therein, for recent applications).

Scattering Amplitude
We consider the scattering amplitude of the process,  reaction are s = (p 1 + p 2 ) 2 , t = (p 1 − p 3 ) 2 , and u = (p 2 − p 3 ) 2 , satisfying the condition s + t + u = 2M 2 . In the physical region, the range of Mandelstam variables reads, (2. 2) The dependence of the scattering amplitude on the kinematic variables can be conveniently parametrised in terms of the dimensionless variables, η and φ, defined as, which, in the physical region satisfy the conditions, The scattering amplitude A of the process can be evaluated in perturbative QCD, and expressed as a power series in the strong coupling α s , as, The LO term A (0) , referred to as Born term, receives contribution from a single tree-level Feynman diagram, see Fig. 1. The colour-summed, un-polarised squared amplitude at LO (summed over the number of colours, summed over the final spins, and averaged over the initial states) has a rather simple expression, with, where N c is the number of colours, and = (4−d)/2, with d being the number of (continuous) space-time dimensions. The higher order contributions A (n) , with n = 1, 2, get contributions and can be organised as combinations of gauge invariant factors, according to the dependence on the number of colours (N c ) and on the flavour structure (i.e., the number of light-and heavy-fermion closed loops, respectively, n l and n h ). In particular, the contributions at oneand two-loop admit the following decomposition [36,37,48], The analytic expressions of the one-loop form factors have been known since long time [1][2][3][4][5][82][83][84].
Regarding the two-loop form factors in the colour decomposition (2.10), contributions from the leading colour (A (2) ), one closed fermionic loop (D h ), and two closed fermionic loops (F h ) are known both numerically as well as analytically [36,39,40]; B (2) and C (2) , instead, are known only numerically [36]. Their analytic evaluation requires the evaluation of non-planar diagrams (that give no contribution to the leading colour term), and they are presented for the first time in this work.
The evaluation of the previously known colour factors, together with the novel calculation of B (2) and C (2) , allows us to obtain, for the first time, the complete analytic expression of the two-loop scattering amplitude for the four-quark scattering in QCD with a massive quark-pair, both as internal and as external states.
The results for the four-quark scattering qq → QQ in QCD, hereby presented, can be considered as the natural extension to a non-Abelian theory of the ones obtained for the four-fermion scattering e + e − → µ + µ − in QED, recently presented in [79]. We observe that the coefficient C (2) , as well as E lh , and F (2) h , can be written as linear combination of (colour stripped) Feynman diagrams that appear also in the Abelian case. The form factors A (2) , B (2) , D (2) l , and D (2) h get contribution from Abelian and non-Abelian (colour stripped) diagrams. We refer the Reader to Appendix A for a detailed discussion on the colour decomposition.
The complete analytic calculation of M (2) , or in other words, the computation of the form factors in decomposition (2.10), is the main result of the present manuscript.

Analytic Evaluation
We begin by considering the bare LO squared amplitude and the bare interference terms, are the coefficients of the series expansion of the bare amplitude A b in the bare strong coupling constant, α b s ≡ g 2 s /(4π). Its expression up-to the second-order corrections reads as, with S ≡ (4πe −γ E ) , and µ being the 't Hooft mass scale. The CDR scheme is adopted throughout the whole computation, hence, internal and external states are, accordingly, regularised in d = 4 − 2 space-time dimensions [85][86][87]. The LO term A , given in Eq. (2.6), is finite in the limit d → 4 ( → 0); whereas the higher order terms require the evaluation of one-and two-loop integrals that may contain UV and IR divergences, parametrised as poles in .

UV Renormalisation
The one-and two-loop amplitudes contain both UV and IR divergences. The UV divergences are removed by renormalising the bare quark fields and the bare mass of the heavy quark in the on-shell scheme,     By employing this, we can express the renormalised amplitude in terms of the bare amplitude as, where Z 2,q and Z 2,Q are the on-shell wave function renormalisation constants for the massless and massive quarks; M is the renormalised mass for the heavy quark in the on-shell scheme.
The renormalised amplitude depends on four renormalisation constants (Z 2,q , Z 2,Q , Z αs , Z M ), which admit a perturbative expansion in the renormalised coupling constant α s , The mass and wave-function renormalisation of the heavy quark is known to three loop accuracy in the on-shell scheme [88][89][90]; the wave-function renormalisation of the light quark, due to the presence of heavy quark, is provided at two loop accuracy in [48]; the strong coupling constant renormalisation is known up-to five-loop accuracy [91][92][93][94][95][96]. Their expressions, up-to the required order, are collected in Appendix B. Upon combining Eqs. (3.3), (3.6), and (3.7), we obtain the UV renormalised amplitude A, given in Eq. (2.5), whose coefficients A (n) can be written in terms of the bare coefficients A (n) b , as, with, The last term in Eq. (3.9b), corresponding to the mass renormalisation counter-term, takes contributions from the diagrams depicted in Fig. 5 and consists of the one-loop diagrams with an insertion of the mass counter-term in the heavy-quark propagators.
With the above definitions, one-and two-loop renormalised interference terms M (n) are obtained as, where, (3.11)

Algebraic decomposition
The generation of the one-and two-loop diagrams contributing to M (1) b and M (2) b , as well as of those needed for the mass-renormalisation, is carried out using FeynArts [52]. By choosing Feynman gauge for the gluon propagator, we identify 10 diagrams at one loop, 184 diagrams at two loops, and 7 counter-term diagrams for the mass renormalisation, respectively, shown in Fig. 2, Figs. 3 and 4, and Fig. 5. Scaleless loop integrals (e.g., oneand two-loop massless tadpoles), and non-planar diagrams that vanish because of colour algebra (see Fig. 6) are neglected. 1 Figure 6: Two-loop diagrams that, upon interference with the Born amplitude, give rise to vanishing contributions, due to colour algebra.
After performing colour, spin and Dirac-γ algebra by means of FeynCalc [53], the interference terms are expressed in terms of n-loop scalar integrals as, where G denotes an n-loop graph interfered with the Born terms, D σ denotes the set of denominators corresponding to the internal lines of G, and N G stands for a polynomial in the scalar products built out of external momenta p i and loop momenta k i , and M 2 . The decomposition of the integrals is automated within the Aida framework [49], where integrands are grouped according to their common set of propagators with respect to the ones of the parent graphs, identified among all the diagrams as the ones with the largest sets of independent denominators. At one-loop, Aida identifies 3 parent graphs, shown in Fig. 7; at two-loop, 31 parent diagrams (22 belonging to four-point topologies and 9 belonging to three-point topologies), shown in Fig. 8, for representative topologies.
The quantities M (n) are simplified within Aida by employing the adaptive integrand decomposition method [50,51] followed by the use of integration-by-parts identities [97][98][99]. The latter are automatically generated for the parent diagrams only, generated by Aida through its interface to the public codes Reduze [54] and Kira [55]. After integrand and integral decompositions, the interference terms M  [102], and independently in [65]. The one-loop counter-term δM (1) is directly computed from the knowledge of the renormalisation constants δZ j and the Born squared amplitude. Differently, the two-loop counter-term δM (2) requires also the decomposition of one-loop integrals, due to both the genuine one-loop amplitude and to the mass renormalisation counter-term, coming from the one-loop diagrams shown in Fig. 5, and, therefore, it admits a decomposition in terms of the basis I (1) .

Results
After inserting the expression of the MIs and adding the bare quantities M to the corresponding counter-terms δM (n) , finally, the renormalised interference terms M (n) are analytically expressed as a Laurent series around = 0, as whose coefficients M (n) k contain GPLs, iteratively defined as [103], The analytical expression of M (1) and M (2) are computed in the non-physical region, s < 0, t < 0, and their analytic continuation to the region of heavy-quark pair production is performed numerically. In particular, M (2) contains 5033 GPLs up-to weight four, whose arguments are written in terms of 18 letters, w i = w i (x, y, z), which depend on the Mandelstam variables through the relations [63,64,100], The numerical evaluation of GPLs, in the physical region (2.4), is performed by adopting the prescription, s → s + iδ , (4.4) by assigning a small positive imaginary part to the squared center-of-mass energy variable, above the pair production threshold. 3 As anticipated in Sec. 2, the analytic evaluation of the one-loop amplitude has been performed long ago by following different approaches [1][2][3][4][5][82][83][84]. On the two-loop side, instead, analytic expressions for the form factors present in the colour decomposition (2.10) is partially known. In particular, the knowledge of these analytic expressions is restricted to leadingcolour and closed fermion-loop terms (A (2) , D (2) l , D (2) h , E (2) l , E (2) h , F (2) l , F (2) lh , F (2) h ) [39,40]. The analytic evaluation of B (2) and C (2) required the evaluation of non-planar diagrams, that were absent from the leading-colour term, and they are presented for the first time in this work.
The independent evaluation of the previously known form factors, together with the novel calculation of B (2) and C (2) , allows us to validate the previously known numerical results [36], and to obtain, for the first time, the complete analytic expression of the two-loop scattering amplitude for the partonic scattering qq → QQ in QCD. Our result is the first example of a complete analytic calculation of a two-loop amplitude in QCD with a massive quark-pair in the internal and as well as external states, including both the leading and sub-leading colour contributions.
A flow chart of the complete computational algorithm implemented in the Aida package is shown in Fig. 9.

IR structure
The structure of IR singularities of the massless and massive gauge theory scattering amplitudes has been studied in [66][67][68][69][70][71][72][73][74][75][76]. The coefficients of the poles in appearing in the renormalised amplitudes M (1) and M (2) agree with the the universal IR structures of the QCD amplitudes, derived from the knowledge of the lower order terms, within SCET [74,[76][77][78]104], poles , (4.5b) where Z IR i (i = 1, 2) are the coefficients of the IR renormalisation factor Z IR encoding the IR divergence. For the process under consideration, involving the production of a massive quark pair, Z IR reads as [78], with, where Γ i = ∂Γ i /∂ ln µ, and Γ i and β i are the coefficients of the perturbative expansion of the anomalous dimensions and of the QCD beta-function, respectively (expressed in the powers of renormalised coupling constant α s ). The anomalous dimension matrix Γ for the qq → tt has been reported in [78].

Finite terms
In Fig. 10, we plot the finite part of one-and two-loop renormalised amplitudes M The plots are obtained by evaluating the analytic formulas at one and two loops with high precision on 10,500 evenly spaced grid points. The numerical evaluation of the GPLs is carried out by HandyG [59] (away from threshold) and Ginac [58] (close to threshold) through their interface to PolyLogTools [57].
The finite term of the analytic expression of the two-loop contribution M 0 , which constitutes the main result of this communication, is found in agreement with the numerical results available in [36]. In particular, the numerical values of the grid attached to the arXiv submission of the latter reference agree with the (higher accuracy) values obtained from the numerical evaluation of our analytic expressions, in the same phase-space points. For completeness, the values of M (2) 0 , numerically evaluated at 1600 phase-space points, are given in the ancillary file qqQQGrid.m, attached to this publication. Our grid is given in the format {φ, η, M 0 }, for the scale choice µ 2 = M 2 , in the same phase-space points chosen in [36].
In Table 1, we showcase the numerical values of the analytic expressions of the individual colour factors, at one-and two-loop. The analytic expressions are evaluated with Ginac at the kinematic point s/M 2 = 5, t/M 2 = −5/4, µ 2 = M 2 (following the prescription given in eq.(4.4), the imaginary term is chosen to have δ = 10 −25 ), which corresponds to the same kinematic point as in the Table 3 of [36] (see also Table 1 of [37]), and our results are in agreement up-to the digits reported in the latter.
Moreover, the analytic expressions of the finite part, as well as of the poles, of A (2) , D h agree with earlier results published in [39,40].    Table 1 of [37].

Conclusion
We completed the analytic evaluation of the scattering amplitude for the process qq → QQ at two loops in QCD, for a massless (q) and a massive (Q) quark type. The contribution of the leading colour diagrams and of those containing fermion loops, whose analytical results were already available in the literature, were independently evaluated and cross checked, and combined with the novel contributions of the sub-leading colour terms, which were evaluated in this work, for the first time.
The un-renormalised interference terms of the one-and two-loop bare amplitudes with the leading-order one were computed in the framework of CDR. The renormalisation of the ultraviolet divergences was carried out by employing the on-shell scheme for the quarks, and the MS-scheme for the strong coupling constants.
The analytic results of the one-and two-loop renormalised contributions, obtained as Laurent series around d = 4 dimensions, respectively, up-to the first order term, and up-to the finite term, were expressed in terms of GPLs and transcendental constants of up-to weight four. The singularity structure of the renormalised results was found to be in compliance with the predicted universal infrared behaviour of QCD amplitudes [76][77][78]. Numerical and partial analytical results of the scattering amplitude already available in the literature [36,37,39,40] agree with the novel analytic expression.
The analytic results of the two-loop scattering amplitude for the top-quark pair production from the light-quark annihilation channel are an essential ingredient to be combined with the ones of the gluon fusion channel, whose analytic knowledge is partially available [41-44, 46, 47], to obtain -hopefully, in a not-so-far future -the full analytic expressions of the scattering amplitudes for the production of a heavy quark-antiquark pair in hadron collisions, at two loops in QCD [36,37].
The results presented for the process qq → QQ in QCD can be considered as an extension to the non-Abelian case of the ones recently obtained for the process e + e − → µ + µ − in QED [79]. The automatic framework which was developed for these calculations is flexible and applicable to other scattering reactions. The computational efforts and the intermediate results for the non-Abelian case, such as diagram generation, integral and integrand decompositions, and evaluation of master integrals, are ingredients that are now available for the study of the elastic scattering processes of one massless and one massive particle/body, which is related by crossing symmetry to the one presented here.
The competences acquired during this work, as well as the building blocks of the calculations, are not limited to applications within Particle Physics, and could be applied, for instance, to investigate processes in General Relativity, like the bending of light caused by a massive astrophysical body, see for instance [105,106], where the massless quark is replaced by a photon, the massive quark is replaced by the world-line of a black-hole, and gluons are replaced by gravitons.
Let us finally remark that, more generally, the presented results constitute a crucial reference for the study of the scattering of particles/bodies with non-vanishing masses, for interactions mediated by self-interacting massless quanta, in the limiting case when one of the body can be treated as massless. Therefore, they can offer additional insights for investigating similarities and differences between fundamental interactions occurring in different physical scenarios.

A Colour Stripped Form Factors
The Feynman diagrammatic approach has been adopted throughout the calculation, and in this Appendix, we provide details on the contribution of the one-and two-loop Feynman diagrams to the form factors present in decompositions (2.9) and (2.10), respectively.
In decomposition (2.9), for the one-loop contribution, we need to deal with 10 nonvanishing Feynman diagrams (see Fig. 2). Two of them contain vacuum polarisation insertions (with a closed heavy-and light-quark loop) contributing to form factors C (1) l and C (1) h . The remaining 8 diagrams may contribute to either A (1) (5 diagrams) or B (1) (4 diagrams). In particular, A (1) gets contribution from purely planar diagrams with and without self-gluon interactions. B (1) , C (1) l , and C (1) h get contribution only from diagrams without self-gluon interactions.