Top quark pair production at NNLO in the quark-antiquark channel

We present the derivation of the NNLO two-parton final state contributions to top pair production in the quark-antiquark channel proportionnal to the leading colour factor Nc2. Together with the three and four-parton NNLO contributions presented in a previous publication, this enables us to complete the phenomenologically most important NNLO corrections to top pair hadro-production in this channel. We derive this two-parton contribution using the massive extension of the NNLO antenna subtraction formalism and implement those corrections in a parton-level event generator providing full kinematical information on all final state particles. In addition, we also derive the heavy quark contributions proportional to Nh. Combining the new leading-colour and heavy quark contributions together with the light quark contributions derived previously, we present NNLO differential distributions for LHC and Tevatron. We also compute the differential top quark forward-backward asymmetry at Tevatron and find that our results are in good agreement with the measurements by the D0 collaboration.


JHEP12(2015)074
1 Introduction The precise study of top quark properties provides a detailed probe of the Standard Model of particle physics and models beyond it. Due to its small lifetime, the top quark decays before it hadronizes and it is the only quark whose production dynamics can be studied without having to account for hadronization effects. Following its discovery twenty years ago [1,2], the top quark has been studied intensively at the Tevatron. Owing to the limited number of top quarks produced, differential studies there suffered from rather large uncertainties [3,4]. The measured top quark forward-backward asymmetry received much attention since it was first found to deviate substantially from Standard Model expectations [5].
With the large number of top quark pairs produced at the LHC, the study of its properties has become precision physics, and the measurement of very precise differential cross sections, attainable. Recently, the ATLAS and CMS collaborations reported measurements of normalised differential observables in tt production in several kinematical variables, such as the transverse momentum, rapidity and invariant mass of the tt system, as well as the top quark transverse momentum and rapidity [6][7][8][9][10]. More recently, both LHC collaborations have measured the charge asymmetry in top-pair production, [11,12].
These measurements will allow for a very detailed and accurate probe of the top quark production mechanism. To reliably interpret data, these very precise measurements must be matched onto equally accurate theoretical predictions, which can be obtained by including contributions up to next-to-next-to leading order (NNLO) in perturbative QCD. At present, NNLO corrections have been included in calculations of the inclusive total tt cross section in [13], for the inclusive and differential Tevatron top quark forward-backward asymmetry in [14], and for differential LHC distributions in [15].
At NNLO, perturbative calculations of collider observables are typically carried out using parton-level event generators. These programs generate events for all parton-level subprocesses relevant to a given final state configuration up to NNLO accuracy and provide full kinematical information on an event-by-event basis. An NNLO event generator for observables with n final state particles or jets involves three main building blocks: the two-loop corrections to the n-parton final state, referred to as double-virtual contributions (dσ VV NNLO ), the one-loop corrections to the (n + 1)-parton final state, called real-virtual contributions (dσ RV NNLO ), and the tree-level (n + 2)-parton double real contribution (dσ RR NNLO ). These three building blocks involve infrared divergences that arise from the exchange or emission of soft and collinear partons and cancel only in their sum. In addition, two mass factorisation counter-terms, dσ MF,1 NNLO and dσ MF,2 NNLO , are needed in the three and two-parton final states contributions respectively in order to cancel infrared divergences originated from initial state collinear radiation.
The combination of subprocesses of different particle multiplicity which are individually infrared divergent is a major challenge in the construction of NNLO parton-level event generators. Employing a subtraction method to regulate these infrared singularities, the NNLO partonic cross section for top pair production in a given partonic channel has the -2 -JHEP12(2015)074 general structure [16] dσ NNLO Two types of subtraction terms are introduced: dσ S NNLO for the 4-parton final state, and dσ VS NNLO for the 3-parton final state. The former approximates the infrared behaviour of the double real contributions dσ RR NNLO in their single and double unresolved limits, whereas the latter reproduces the single unresolved behaviour of the mixed real-virtual contributions dσ RV NNLO .
In the context of the present computation, we further decompose the double real subtraction term dσ S NNLO as in [17], which allows us to rearrange the different terms in eq. (1.1) into the more convenient form The NNLO contributions to the top pair production in the quark-antiquark channel can be decomposed into colour factors as 5) with N c being the number of colours, N l the number of light quark flavours, and N h the number of heavy quark flavours. All coefficients multiplying the different colour factors in eq. (1.5) are gauge invariant and can be computed independently. The goals of this paper are twofold. The first aim is the derivation of the two-parton contribution given in eq. (1.2) for the colour factor N 2 c within the theoretical framework of antenna subtraction with massive fermions. This framework has been established in the massless case in [16,[18][19][20][21][22][23] and extended to the massive case in [17,[24][25][26][27][28][29]. We

JHEP12(2015)074
shall not repeat it here. Together with the three and four-parton contributions derived in [17] the two-parton contribution yields the coefficient A in eq. (1.5) and enables us to obtain theoretical predictions for the most significant NNLO contribution to top pair hadro-production in the quark-antiquark channel. In addition, we derive the heavy flavour contributions labeled in eq. (1.5) as F h and G h , presenting a detailed discussion of the ultraviolet renormalisation of the loop amplitudes involved.
Our second goal is to present phenomenological results for Tevatron and LHC energies, including the leading-colour and heavy quark contributions derived in this paper, as well as the light quark terms, which we completed in previous papers [28,29].
We employ analytic expressions for all matrix elements, subtraction terms and integrated subtraction terms that arise in our calculation. For the real-virtual contributions, we use OpenLoops [30] in combination with CutTools [31] and for the double virtual contributions we use the analytic two-loop matrix-elements of [32,33]. 1 As a result, at the real-virtual and virtual-virtual levels, the explicit infrared poles are cancelled analytically. Furthermore, the numerical evaluation of the finite remainders at each level included in the parton-level generator is substantially faster than when these terms are evaluated numerically.
The paper is organized as follows: in section 2 we present the derivation of the new integrated massive initial-final antennae required in the double virtual counter-term dσ U NNLO for the leading-colour contributions to top pair production in the quark-antiquark channel. Section 3 contains the explicit derivation of the counter-term dσ U NNLO expressed in terms of massive integrated dipoles. Employing those dipoles, we find that the virtual-virtual subtraction term has a similar structure to the one found in the massless case in the context of di-jet production at NNLO [18]. Section 4 is dedicated to the heavy quark contributions proportional to N h N c and N h /N c . Section 5 contains our numerical results: various differential distributions are presented and the phenomenological impact of this computation is discussed. Finally section 6 contains our conclusions. We enclose three appendices containing the master integrals required to compute the integrated antennae presented in section 2, as well as the phase space parametrisations needed to calculate these integrals.
2 Integrated initial-final massive antennae

General features
One of the aims of this paper is the computation of the NNLO two-parton contribution proportional to N 2 c for top pair production in the quark-antiquark channel. Within the antenna subtraction method, this involves the construction of the double virtual counterterm dσ U qq,NNLO,N 2 c which renders the two-parton contribution finite and suitable to be implemented in an NNLO parton-level generator.

JHEP12(2015)074
As shown in eq. (1.4), this virtual-virtual counter-term contains mass factorisation terms together with integrated double real and real-virtual subtraction terms. In order to compute these, the integrated forms of the antenna functions present in the corresponding real-virtual and double real subtraction terms are needed. The integrated antennae are obtained by integrating the antenna functions over the corresponding antenna phase space inclusively. This integration is carried out analytically, and turns the implicit soft and collinear singularities of the antenna functions into explicit poles in the dimensional regularisation parameter , which are then cancelled against the explicit poles of the double virtual matrix elements and mass factorisation counter-terms.
In the context of this paper, we need the integrated forms of two massive initial-final A-type antennae which are presented here for the first time. These antennae involve a massless initial state quark and a massive quark in the final state as radiators, and one or two final state gluons which are unresolved, i.e. soft and/or collinear. More concretely, the calculation that we present in this paper requires the integrated forms of the treelevel four-parton antenna A 0 4 (1 Q , 3 g , 4 g ,2 q ) and the one-loop leading-colour three-parton antenna A 1,lc 3 (1 Q , 3 g ,2 q ). The tree-level four-parton antenna appears in the double real subtraction term (dσ S,2 qq,NNLO,N 2 c ) given in [17]. It is needed in order to capture the infrared behaviour of the leading-colour double real contributions associated to the partonic process qq → ttgg when the two final state gluons are unresolved and colour-connected. The oneloop antenna appears in the real-virtual subtraction term dσ VS,a qq,NNLO,N 2 c given in [17]. It is required to capture the infrared behaviour of the leading-colour part of one-loop matrix element associated to the partonic process qq → ttg.
The general definition of the integrated forms of three and four-parton massive antennae in all three configurations (final-final, initial-final, initial-initial), as well as the corresponding phase space factorisations and mappings, have been presented in [24,27,28]. Below we will only briefly recall the definition of initial-final integrated massive antennae, as they are the only ones needed in the calculation presented in this paper.
Initial-final three-parton tree-level and one-loop antennae like A l 3 (1 Q , 3 g ,2 q ) (l = 0, 1) are generically denoted as X l i,jk , with parton i in the initial state, and j, k in the final state. The kinematics associated to these initial-final antenna functions with one massive radiator are p i + q → p j + p k , with q 2 = −Q 2 < 0, p 2 i = p 2 j = 0 and p 2 k = m 2 Q , and the integrated antennae are defined as [28] dΦ 2 is the corresponding 2 → 2 phase space, and 3) The initial-final massive antenna phase space denoted by dΦ X i,jk is given by

JHEP12(2015)074
In addition to the phase space integration, one-loop antennae must be also integrated over the loop momentum. Initial-final four-parton antennae like A 0 4 (1 Q , 3 g , 4 g ,2 q ) are generically denoted as X 0 i,jkl , with parton i in the initial state and j, k, l in the final state. The associated kinematics is p i + q → p j + p k + p l , with q 2 = −Q 2 , p 2 i = p 2 j = p 2 k = 0 and p 2 l = m 2 Q . The integrated forms of these four-parton initial-final antennae are obtained as 5) and the corresponding antenna phase space is given by After the inclusive phase space integration, the three and four-parton integrated antennae are functions of Q 2 , p i · q and m 2 Q . We change the dependence on the latter two variables for and We also change the dependence on Q 2 for sī ( jk) and sī ( jkl) in three and four-parton integrated antennae respectively, where the bars and tildes are related to the initial-final phase space mappings that can be found, for example, in appendix B of [28]. The invariants are defined as s ij = 2p i · p j with the crossings explicitly performed on the momenta p i , p j , i.e. p 0 i , p 0 j > 0. In integrated subtraction terms, the mapped final state momenta p jk and p jkl are relabeled in accordance to the labeling of momenta in the lower multiplicity final state, in such a way that the integrated antennae end up depending on invariants sī j . This is the notation that we shall follow throughout this paper.
In all cases, the inclusive phase space integration is carried out following the standard technique of reduction to master integrals using integration-by-parts identities (IBP) [36,37]. The computation of the master integrals is then carried out either directly to all-orders or iteratively, order by order in , using differential equations techniques [38].
All the master integrals found contain multiplicative factors of the form (1 − x i ) −n which regulate soft endpoint singularities and should be kept unexpanded. All other terms can be safely expanded in . As we explain in detail in appendix C, in some of the master integrals found in the reduction of the one-loop antenna of interest here, namely A 1,lc 3 (1 Q , 3 g ,2 q ), more than one -power of (1 − x i ) is encountered. Therefore, the master integrals which we denote collectively by I α (x i , x 0 , ) can be most generally expressed as: The integers m and n are specific to each master integral; the functions R α (x i , x 0 , ) are regular as x i → 1 and can be calculated as Laurent series in .

JHEP12(2015)074
The integrated antennae collectively denoted by X (x i , x 0 , ) are linear combinations of master integrals with coefficients containing poles in as well as in (1 − x i ). After the masters have been inserted into the integrated antennae, these take the form where R X (x i , x 0 , ) is a regular function as x i → 1. The Laurent expansion of the singular factor (1 − x i ) −1−n is done in the form of distributions: It is worth noting that in the functions R α (x i , x 0 , ) and R X (x i , x 0 , ), which are regular as x i → 1, the massless limit x 0 → 1 cannot in general be safely taken. This is due to the presence of terms of the form log k (1 − x 0 ) = log k (m 2 Q /(Q 2 + m 2 Q )), which are expected and correspond to quasi-collinear limits of the antenna functions [24,27,39].
In order to perform the Laurent expansion of the integrated antennae, we distinguish two regions: a hard region where x i = 1 and a soft region where x i = 1. The highest order in needed in the expansion of each master integral is determined by the and x i -dependent coefficient that multiplies the integral in the integrated antenna. In the soft region, since the expansion in distributions generates an additional 1/ factor, the functions R (n) α are required to one order higher in than in the hard region. Further details concerning the integration of A 0 4 (1 Q , 3 g , 4 g ,2 q ) and A 1,lc 3 (1 Q , 3 g ,2 q ) will be presented below.

2.2
The integrated tree-level antennae A 0 q,Qg and A 0 q,Qgg The computation of dσ U qq,NNLO,N 2 c requires the knowledge of the integrated tree-level three and four-parton massive initial-final antennae A 0 q,Qg and A 0 q,Qgg . The former has been derived in [24], and only its pole part will be recalled below. The integrated flavourviolating four-parton antenna A 0 q,Qgg is new and will be presented here for the first time.

The integrated three-parton antenna
The pole part of the integrated antenna A 0 q,Qg is given by where I Qq is a colour-order infrared singularity operator which has been presented in [24] and reads -7 -

JHEP12(2015)074
and Γ (1) qq (x) is a colour-ordered splitting kernel given by As can be seen from the equations above, A 0 q,Qg has its deepest pole at 1/ 2 . In the virtual-virtual subtraction term that will be presented in section 3, we will need products and convolutions of two of these integrated antennae, and we will therefore require its Laurent expansion up to order 2 . The resulting expression is too long to be presented here. It is included in an ancillary Mathematica file attached to the arXiv submission of this paper.
The integrated initial-final antennae depend explicitly on one momentum fraction x i carried by the initial state parton. In order to combine these integrated forms with other integrated subtraction terms depending on both momentum fractions x 1 and x 2 carried by the incoming quark-antiquark pair, it is useful to make all antennae explicitly depend on both x 1 and x 2 . We achieve this by introducing delta functions in the following way: The same prescription will be applied to A 0 q,Qgg and A 1,lc q,Qg .

The integrated four-parton antenna
To compute the integrated initial-final massive antenna A 0 q,Qgg we start by expressing all phase space integrals as cuts of two-loop four-point functions with two off-shell legs in forward scattering kinematics [40], and reduce these four-point functions to master integrals with LiteRed [41].
As mentioned above, we consider the following DIS-like kinematics The denominators that we find are: where D 1 , D 2 and D 3 are the cut propagators. In the reduction procedure we impose momentum conservation and mass-shell conditions for the external legs, and we discard all (b) I [4] .
(e) I [4,7,8] . integrals that do not contain the propagators D 1 , D 2 or D 3 with power -1. We find the six master integrals depicted in figure 1. Their explicit expressions, together with a parametrisation of the phase space measure are presented in appendices A and B respectively.
The integrals I [0] , I [8] and I [4] are known [24,25]. We computed the remaining master integrals, namely, I [4,8] , I [4,5,8] and I [4,7,8] , using differential equation techniques, and we present them here for the first time. Although the differential equations for these three integrals are coupled, they can be decoupled order by order in , which enabled us to solve them iteratively. We fixed two of the three sets of undetermined integration constants by demanding that the integrals be regular as x → 1 after factoring out an appropriate power of (1 − x). In order to fix the remaining set of integration constants we calculated via direct evaluation the soft limit of I [4,7,8] , and employed it as a boundary condition. The details of the computation of this soft limit can be found in appendix B.
The highest order in with which each master integral contributes up to the finite part of A 0 q,Qgg is specified in table 1. As mentioned before, in the soft region (x = 1), they must be expanded to one order higher than in the hard region (x = 1).
After expanding in distributions we find that the integrated four-parton massive antenna A 0 q,Qgg has a deepest pole at O(1/ 4 ). This is expected since at the unintegrated level A 0 4 (1 Q , 3 g , 4 g ,2 q ) contains double soft and triple collinear singular limits, which are employed to reproduce the singular behaviour of the double real contributions associated to the partonic process qq → ttgg. The integrated antenna can be written in terms of harmonic polylogarithms (HPLs) with arguments x i or x 0 and generalised harmonic polylogarithms (GPLs) of argument x 0 and weights involving 1/x i . GPLs and HPLs appear with up to trascendentality three and four respectively. Given its length, the full expression of this antenna up to O( 0 ) is given in the ancillary file attached to the arXiv submission -9 -

JHEP12(2015)074
Integral Needed in the hard region to Deepest pole of this paper. The deepest poles read The leading-colour one-loop antenna employed in [17] for the construction of a real-virtual subtraction term for partonic process qq → ttg at one-loop is given by where the subscript [lc] indicates that only the leading-colour primitive amplitude (proportional to N c ) is kept. As mentioned above, the DIS-like kinematics associated with this antenna function is The corresponding phase space parametrization has been derived in [28] in the context of the integration of tree-level initial-final massive antennae. It is repeated for completeness in appendix C.
The phase space integration of the second term in eq. (2.21) is trivial, since the ratio |M 2 (( 13) Q ,2 q )| 2 /|M 0 2 (( 13) Q ,2 q )| 2 only depends on Q 2 , and it can be therefore pulled out of the phase space integral. Thus, only A 0 3 (1 Q , 3 g ,2 q ) must be integrated over the phase space, and this integral is known.
(g) I [3,4,6,7] . We must therefore focus on the first term in eq. (2.21), for which we follow the same procedure as the one described above for A 0 q,Qgg . The following denominators must be considered: where D 1 and D 2 are the cut propagators, D 3 to D 6 are the loop propagators, and D 7 = −s 23 . We find the master integrals depicted in figure 2.
In addition to these new master integrals, the ultraviolet renormalised integrated antenna A 1,lc q,Qg also contains the known inclusive phase space integral, denoted by

JHEP12(2015)074
Integral Needed in the hard region to Deepest pole was computed in [24]. This integral is brought about by the phase space integration of the mass and strong coupling renormalisation counter-terms.
The evaluation of the loop integrals I [3] , I [3,5] , I [3,6] and I [4,6] is rather straightforward. It can be done to all orders in by integrating the underlying loop integrals over the phase space. The remaining three master integrals, namely I [3,4,6] , I [3,4,5,6] and I [3,4,6,7] , cannot be computed in this way. Instead, we calculated them using differential equations. Unlike the case of the master integrals needed for A 0 q,Qgg , the integration constants cannot be in this case fixed by imposing regularity conditions on the integrals. This impossibility is due to the fact that each integral contains several different powers of (1 − x). Therefore, in order to determine the integration constants, we employed in this case independently computed soft limits as boundary conditions. For I [3,4,5,6] we used the soft limit of the underlying one-loop box given in [42], whereas for I [3,4,6] and I [3,4,6,7] we calculated the soft limit of the one-loop triangle to all orders in using a Mellin-Barnes representation. More details on these derivations are given in appendix C.
Like the integrated four-parton antenna presented above, A 1,lc q,Qg has its deepest pole at O(1/ 4 ). It can also be written in terms of HPLs and GPLs of trascendentality four and three respectively. The HPLs have arguments x i or x 0 , and the GPLs have argument x 0 and 1/x i in the weights. The complete expression for A 1,lc q,Qg is enclosed in the ancillary file included in the arXiv submission of this paper. Its deepest poles are given by 3 Top quark pair production at NNLO Following the decomposition of the NNLO cross section for top pair production in two, three and four-parton final state contributions as shown in eqs. (1.2) and (1.4), in this section we present the two-parton final state for the coefficient A in eq. (1.5). In particular, we shall construct the virtual-virtual subtraction term dσ U qq,NNLO,N 2 c as a combination of integrated real-real and real-virtual subtraction terms, and mass factorisation counter-terms. We will show that these ingredients can be arranged in such a way that dσ U qq,NNLO,N 2 c can be expressed in terms of so-called integrated dipoles, which facilitate the cancellation of explicit infrared poles in the two-parton final state contribution.
For massless jet observables computed at NNLO with antenna subtraction, the general formalism allowing to write the real-virtual and virtual-virtual counter-terms using integrated dipoles was presented in [18]. There it was shown that massless integrated dipoles, denoted as J ( ) 2 with = 1, 2, emerge naturally when rearranging the mass-factorisation counter-terms and integrated double real and real-virtual subtraction terms into three different pieces that are free of collinear initial state singularities. It was furthermore observed that the colour-ordered J ( ) 2 's are related to Catani's one and two-loop infrared singularity operators I ( ) ij ( ) [43]. As we will see below, the virtual-virtual subtraction term dσ U qq,NNLO,N 2 c for heavy quark pair production is similar in structure to the massless case. For the purpose of making this similarity manifest, we define a new type of massive integrated dipoles J ( ) 2 . In particular, the dipole J (1) 2 that we find in the present calculation of this leading colour contribution is related to the colour-ordered infrared singularity operator I (1) Qq , defined in [17,24] and recalled in eq. (2.14).

NNLO contributions to top quark pair production in the qq channel
In order to set up our conventions for normalisation factors and matrix elements, we start by recalling the leading order partonic cross section for heavy quark pair production. It can be written as where, dΦ 2 (p 3 , p 4 ; p 1 , p 2 ) is the 2 → 2 partonic phase space, and J (2) 2 (p 3 , p 4 ) is a so-called measurement function, which ensures that a pair of final state massive quarks of momenta p 3 and p 4 are observed. M 0 4 (. . .) is a colour-ordered and coupling-stripped tree-level amplitude. It is related to the full amplitude through the (trivial) colour decomposition

JHEP12(2015)074
The normalisation factor is where s is the energy squared in the hadronic center-of-mass frame. Included in N qq LO are the flux factor, as well as the sums and averages over colour and spin. The constant C( ) has been defined in eq. (2.3), whileC( ) is given bȳ providing the useful relation The leading-colour double real corrections to heavy quark pair production are due to the tree-level partonic process qq → QQgg. They read where dΦ 4 is the 2 → 4 phase space, and |M 0 | 2 is the square of the full coupling-stripped tree-level amplitude normalised to (N 2 c − 1). This factor is included in the overall normalisation N RR,qq NNLO , which is given by The subscript N 2 c on the matrix element squared indicates that only the terms proportional to this colour factor are to be kept.
Given the fact that the measurement function J (4) 2 allows the final state gluons with momentum p 5 and p 6 to become unresolved, i.e. soft or collinear, eq. (3.6) contains infrared divergences. These divergences are implicit, in the sense that they only become explicit as poles in after the integration over the phase space is performed. The antenna subtraction term that regulates these infrared divergences has been derived in [17].
The mixed real-virtual contributions are given by the phase space integral of the interfered one-loop and tree-level amplitudes for the 2 → 3 process qq → QQg. They read with the normalisation factor N RV,qq NNLO given by The real-virtual contributions in eq. (3.8) contain explicit ultraviolet and infrared divergences as well as implicit infrared ones. The explicit poles in originate from the loop integration in M 1 q 1q2 →t 3t4 g 5 , whereas the implicit singularities are due to the phase space integration over regions where the matrix elements diverge: the soft limit p 5 → 0 and the collinear limits p 1 ||p 5 , p 2 ||p 5 . While the ultraviolet poles are cancelled upon renormalisation, we employ a subtraction term dσ T qq,NNLO,N 2 c to deal with the infrared ones. This subtraction term has the twofold purpose of canceling the explicit poles, whose structure is well known [44], while simultaneously regularising the phase space integrand in the soft and collinear limits. Its explicit construction was derived in [17].
Finally, the double virtual contributions have not been derived before and are the object of study in the next section. They involve the interference of a two-loop 2 → 2 matrix element with its tree-level counterpart, and a one-loop amplitude squared. They can be written as, where the normalisation factor N VV,qq NNLO is given by After ultraviolet renormalisation, these double virtual contributions contain explicit infrared poles up to order four, which originate from the loop integration. There are no implicit poles, as the measurement function J (2) 2 does not allow any final state particle to be unresolved. As we shall see below, all explicit infrared poles in eq. (3.10) are captured and cancelled by those in the subtraction term dσ U qq,NNLO,N 2 c . For the two-loop matrix element in eq. (3.10) we employ the analytic results of [33]. The square of the one-loop amplitude was calculated in [45]. However, we independently computed our own expression, written in terms of GPLs. For the numerical evaluation of the HPLs in the double virtual amplitudes we use Chaplin [46], and for the genuine GPLs, the GiNaC implementation of [47]. We derived an expansion of both terms in dσ VV qq,NNLO,N 2 c at the production threshold (ŝ ∼ 4m 2 Q ) in order to circumvent the numerical instabilities that occur when evaluating the GPLs in this limit.
Before proceeding to the construction of the double virtual counter-term, it is worth commenting on the purpose of the x 1 and x 2 dependence in eqs. (3.8) and (3.10). This seemingly trivial dependance introduced in dσ RV In general, the virtual-virtual subtraction term denoted as dσ U NNLO is built from mass factorisation counter-terms, as well as integrated double real and real-virtual subtraction terms. In the present case, it is given by The unintegrated forms of all these terms have been given in [17]. q,Qgg and the one-loop antenna A 1,lc q,Qg . In general, an additional integrated real-virtual subtraction term, related to a double real subtraction term capturing almost-colour-unconnected unresolved limits, and denoted as 1 dσ VS,c NNLO , should be present in eq. (3.12). Its absence in the present case is due to the fact that at the real-real level the partonic process contains only two gluons.
Following [18] we decompose our virtual-virtual counter-term as where, by construction, each term will be free of explicit initial state collinear poles. In order to achieve this, some rearrangements in the integrated real-virtual subtraction terms (3.14) In the first term, which will be included in dσ U,a qq,NNLO,N 2 c , we group the contributions of the form X 0 3 |M 1 4 | 2 , whereas in the second term, which will be a part of dσ U,c qq,NNLO,N 2 c , we put the subtraction terms of the form X 1 3 |M 0 4 | 2 . The integrated form of the ultraviolet-scale-compensating subtraction term given in eq. (8.22) of [17] reads where b 0 is the N c coefficient of the QCD beta function at one loop, i.e. b 0 = 11/6. The square brackets in the second line may be expanded to produce two terms: the terms with the pre-factor (sī j /µ 2 ) − will be part of dσ U,c qq,NNLO,N 2 c , while those terms originated from the factor −1 in the square brackets with be part of dσ U,a qq,NNLO,N 2 c . In order to render each of the three pieces of the virtual-virtual subtraction term in eq. (3.13) free of initial state collinear singularities, we must split the mass factorisation counter-term dσ MF,2 qq,NNLO,N 2 c into three parts. This construction is detailed below. In general, for a given partonic process initiated by two partons i and j with momenta p 1 and p 2 , the mass factorisation counter-term dσ MF,2 ij,NNLO to be included at NNLO at the virtual-virtual level reads [18,22,23], where dσ V kl,NLO , dσ MF kl,NLO and 1 dσ S kl,NLO are the NLO virtual contributions, mass factorisation counter-terms and integrated subtraction terms respectively, and the splitting kernels Γ (1) ij,kl and Γ (2) ij,kl are given by To simplify the construction of the double virtual subtraction term dσ U ij,NNLO , Γ ij,kl can be written as suggested in [18] as ab;kl (x 1 , x 2 ), (3.19) such that, Γ ij;kl (x 1 , The reduced kernel Γ ij (z) is related to the usual Altarelli-Parisi splitting functions [48] as where the individual contributions are given by: ij;kl dσ kl,LO .
where dσ V qq,NLO,Nc is the leading-colour NLO virtual cross section, and with Γ Using the expressions for dσ MF qq,NLO,Nc and 1 dσ S qq,NLO,Nc given in [17], we get The convolution of two functions f (x 1 , x 2 ) and g(y 1 , y 2 ) is defined as in such a way that the convolutions appearing in eq. (3.29) take the following forms The convolution of the form Γ qq ⊗ A 0 q,Qg ( , s) has not been computed before. It is needed up to finite order and reads The complete expression can be found in the ancillary Mathematica file attached to the arXiv submission of this paper. Finally, the mass factorisation counter-term denoted as dσ MF,2,c with the leading-colour two-loop kernel Γ (2) qq;qq (x 1 , x 2 ) defined as in eq. (3.20) and Γ (2) qq (x) = As usual, p qq is given by  (3 Q , 4Q,2q,1 q ) and its tree-level counterpart. The content of the square bracket is manifestly free of initial state collinear divergences as is the whole virtual-virtual subtraction term dσ U,a qq,NNLO,N 2 c .
In analogy to the massless case [18], and using a similar notation, we can reformulate eq. (3.39) in terms of integrated dipoles. Indeed, we can define the integrated massive NLO dipoles as: Using the expression in eq. (2.13) for the pole part of A 0 q,Qg one can immediately see that these massive integrated dipoles of the form J 1 2 are directly related to the massive colour-ordered infrared singularity operator I 2 (p 3 , p 4 ). (3.42) In general, for an n-jet hadron collider observable computed at NNLO, dσ U,a NNLO contains integrated dipoles of three types: final-final, initial-final and initial initial. In our case, for the leading colour contributions, the initial-final dipoles defined in eqs.  is not free of explicit initial state collinear poles. In order to remedy this, we define can be written as a convolution of integrated massive dipoles, with no infrared collinear singularities left: It is worth noting that, in general, the terms that we added to dσ U,b qq,qq (x 1 , x 2 ) |M 0 4 (3 Q , 4Q,2q,1 q )| 2 J

JHEP12(2015)074
The expression in the above equation is not free of initial state collinear poles. It becomes so, however, once we add back the contributions subtracted from dσ U,b qq,NNLO,N 2 c in order to render it free from initial state collinear singularities. Indeed we find that has a similar structure than the one found for dσ U,c NNLO in the massless case, and allows for the definition of new initial-final massive integrated dipoles: In terms of these massive integrated dipoles we have As in the massless case, one can see that the integrated dipoles present here are related to integrated subtraction terms which involve genuine NNLO objects. This feature can be regarded as an additional check on the correctness of the construction of the subtraction terms at real-real, real-virtual and virtual-virtual levels in this NNLO computation.

Explicit pole cancellation
With the explicit expressions included above one can show analytically that each of the virtual-virtual subtraction terms dσ U,a This is an ultimate check that demonstrates that we have correctly implemented the subtractions terms at real-real, real-virtual and virtual-virtual levels for the leading-colour contributions to qq → QQ + X at NNLO.

Heavy quark contributions
In this section we present the heavy quark contributions to the top pair production cross section in the quark-antiquark channel computed at O(α 4 s ). They are proportional to the colour factors N h N c (leading-colour) and N h /N c (subleading-colour), and correspond to the coefficients F h and G h in eq. (1.5).

JHEP12(2015)074
Both F h and G h receive contributions in the two and three-parton final states, and require, in both final states, subtraction terms to cancel explicit and implicit infrared divergences. For all loop matrix elements involved we perform the ultraviolet renormalisation in the so-called decoupling scheme, with the gluon wave functions and the heavy quark mass and wave functions renormalised on-shell, and the strong coupling α s renormalised in the MS scheme with N l active flavours. The two-loop amplitudes employed in our calculation [33], were originally computed with the strong coupling renormalized in the MS scheme with N F = N h + N l active flavours, and therefore required a conversion to the decoupling scheme.
In this section we will explicitly construct the two and three-parton contributions to the N h N c and N h /N c colour factors, and we will show how the amplitudes computed in [33] in the full-flavour scheme can be converted to the decoupling scheme via a finite renormalization of α s .

Real-virtual contributions
In the antenna subtraction framework, the three-parton heavy-quark contribution to top pair production at NNLO in the qq channel is given by Following the colour decomposition of the real-virtual matrix-element as presented in [29] dσ RV qq,NNLO,N h can be written as where we have used the following shorthand notation: (3 Q , 4Q; ;2q, 5 g ,1 q ), (4.4) where the gluon is photon-like.

JHEP12(2015)074
After UV renormalisation in the decoupling scheme dσ RV qq,NNLO,N h is free of explicit -poles. It contains, however, infrared implicit poles: the phase space integral is divergent due to integration over the soft limit p 5 → 0 and the collinear limits p 1 ||p 5 and p 2 ||p 5 , where the matrix elements are singular. The subtraction term dσ T qq,NNLO,N h in eq. (4.1) regularises these implicit divergences of the phase space integral. We shall present it below.

Real virtual subtraction term dσ T qq,NNLO,N h
Both the leading-colour (N h N c ) and subleading-colour (N h /N c ) real-virtual subtraction terms only contain terms commonly known in the antenna subtraction literature as dσ VS,a (1) . They are built as products of tree-level antennae and reduced one-loop matrix elements squared. Unlike in the most general case, there are no integrated double real subtraction terms, no mass factorisation counter terms dσ MF,1 , nor subtraction terms involving one-loop antennae. The absence of one-loop antennae in dσ T qq,NNLO,N h is due to the fact that the N h coefficient of the antennae A 1 3 , which would be in principle be required, vanishes in the decoupling scheme. It is also worth mentioning that in the decoupling scheme the reduced one-loop matrix element needed in dσ T qq,NNLO,N h , i.e. the N h part of the one-loop amplitude for the process qq → tt, is finite.
Explicitly, we find, Since, as we mentioned above, |M
Being a pure subtraction term, dσ T qq,NNLO,N h must be added back in integrated form to the cross section at the two-parton level. The required A-type integrated massive and massless antennae are all known [24,27].

qq,NNLO,N h
In this section we present the heavy quark NNLO two-parton contributions to qq → tt in the decoupling scheme, which are given by We shall focus in particular on the construction of the double-virtual counterterm dσ U qq,NNLO,N h , split into a leading-colour (N h N c ) and a subleading-colour (N h /N c) part. The contributions from the double-virtual matrix elements can be written as with N VV,qq NNLO given in eq. (3.11), and with the abbreviation For the N h part of the two-loop matrix element in eq. (4.8) we employ the analytic results of [32] converted to the decoupling scheme via a finite renormalisation of α s that will be described below. The "one-loop squared" term has been computed analytically in [45]. We re-derived it ourselves, also analytically, and use our own result in our event generator. We further compared these results with those provided by Roberto Bonciani in both renormalisation schemes and found full agreement. After UV-renormalisation, dσ VV qq,NNLO,N h contains explicit infrared poles, that are cancelled by the virtual-virtual subtraction term dσ U qq,NNLO,N h derived below.

Virtual-virtual subtraction term dσ U qq,NNLO,N h
Following the decomposition of the mass factorisation counter term dσ MF,2 and dσ U explained in section 3, we find that dσ U qq,NNLO,N h only contains terms denoted there as dσ U,a , involving only the mass factorisation counterterm dσ MF,2,a qq,N N LO,N h and the integrated form of the real-virtual counter term derived in eq. (4.5).
The mass factorisation counterterm is given by

JHEP12(2015)074
The leading-colour part is given by The terms present in the square bracket are manifestly free of initial-state collinear divergences as is the whole double virtual subtraction term dσ U qq,NNLO,N h Nc . Using the same integrated massive NLO dipoles as in section 3.4 for the leading-colour double virtual counter-term contributions, we can rewrite eq. (4.11) as: The subleading-colour heavy quark counter term dσ U qq,NNLO,N h /Nc reads, This subtraction term can also be written in terms of integrated dipoles. In addition to the initial-final massive dipoles defined in eqs. as well as a massless initial-initial dipole defined in [18] and given by In terms of these three type of integrated dipoles, the virtual-virtual subtraction term dσ U qq,NNLO,N h /Nc takes the following form:

JHEP12(2015)074
The pole parts of each of the integrated dipoles in eqs. (4.12) and (4.16) is given by a specific colour-ordered infrared singularity operator I (1) ij given in [17]. In terms of those operators we find that the pole part of the entire double-virtual heavy-quark counter term reads: is finite. It is furthermore proportional to its tree-level counterpart: with the factor R h in the equation above given by and with β = 1 − 4m 2 Q /s 12 . Using this expression for the one-loop amplitude, we can also express the pole structure of the "one-loop squared" contributions in dσ VV qq,NNLO,N h given in eq. (4.7) in terms of infrared singularity operators. For the leading-colour part we find Taking the pole part of the integrated subtraction terms given above (where we have used momentum conservation to set s 24 = s 13 and s 14 = s 23 ) enabled us to verify the pole structure of the two-loop amplitudes in the decoupling scheme. We conclude this section by showing how the one and two-loop amplitudes computed in [33] in the full-flavour scheme can be converted to the decoupling scheme with N l active flavours. This conversion can be achieved with the well known finite renormalisation of α s (see e.g. [49,50]) We can use eq. (4.22) to relate an amplitude M computed in the full theory to its counterpart in the decoupling scheme, denoted as M. We have, such that for the one and two-loop amplitudes we find respectively,

Numerical results
Together with the double real and real-virtual leading-colour contributions derived in [17], the two-parton channel presented in section 3 completes the calculation of the coefficient A in eq. (1.5). We implemented this coefficient along with F l , F h , G l , G h , H l , H h and H lh , in a Monte Carlo parton-level event generator based on the set up of eq. (1.2). While the inclusion of the coefficients H X is rather straightforward, as they are all infrared finite and only receive contributions from double virtual matrix elements, the coefficients F X and G X are non-trivial. F l and G l were presented in [29], and F h and G h were derived in this publication in section 4.
Using an adaptation of the phase space generator employed in the NNLO di-jet calculation of [51], our program, written in Fortran, provides full kinematical information on an event-by-event basis, thus allowing for the evaluation of differential distributions for top pair production including exact NNLO results for the quark-antiquark channel. Naturally, the calculation of total cross sections is also possible.
In this section, we present numerical results for LHC and Tevatron. For LHC, with √ s = 8 TeV, we show differential distributions in the transverse momentum of the top  quark p t T , the top quark rapidity y t , the invariant mass of the tt system m tt and its rapidity y tt . For Tevatron we show differential cross sections in p t T , y t and m tt , as well as in the absolute value of the top quark rapidity |y t |. In a separate sub-section we present our results for the forward-backward asymmetry A F B , differential in the rapidity difference |∆y tt | = |y t − yt| as well as in m tt and p tt T . All the results presented below include all partonic channels at LO and NLO. At NNLO, the contributions to the quark-antiquark channel from the colour factors N 2 We use the pole mass of the top quark m t = 173.3 GeV, and the PDF sets MSTW2008 68cl. The factorisation and renormalisation scales are set equal to the top quark mass µ R = µ F = m t throughout. Where provided, scale variations correspond to the range m t /2 ≤ m t ≤ 2m t .
The total hadronic cross section in the quark-antiquark channel is presented in table 3 for Tevatron, and LHC with center-of-mass energies 7, 8 and 14 TeV. From the comparison with the full NNLO calculation of [52] it can be seen that, due to the omission of subleadingcolour pieces at O(α 4 s ), our results consistently overestimate the NNLO correction to the total cross section, indicating that the combined contribution of the coefficients C and E in eq. (1.5) is negative and non-negligible. It should be noted that the size of these subleading-colour terms decreases with the hadronic center-of-mass energy, as can be seen from the fact that the discrepancies between or results and those of [52] at NNLO are 5% for Tevatron, 3.25% for LHC 7, 3% for LHC 8 and 2.25% for LHC 14.

Differential distributions for LHC
In figure 3 we present differential distributions for top pair production in pp collisions with √ s = 8 TeV at LO, NLO and NNLO, together with the corresponding k-factors. As expected, the impact of the NNLO corrections in the qq channel on LHC cross sections is in general mild, given the dominance of the gluon-gluon initiated process in this scenario. From the ratios NNLO/NLO in the lower panels of figures (a) and (c), it can be seen that the NNLO corrections decrease the p t T and m tt distributions over the entire spectra considered. The decreases range between 3% and 7%, being in both cases more pronounced in the tails of the distributions.

JHEP12(2015)074
In the distributions in y t and y tt , given in figures (b) and (d), it can be seen that in comparison with the full NLO result, NNLO corrections in the quark-antiquark channel shift the cross sections downwards in the central region by 5%. In the very forward and backwards ends of the spectra the impact of these corrections is more substantial, causing an upwards shift of 12% in the y t distribution, and over 50% in the y tt distribution.
All plots in figure 3 show a slight reduction in the scale uncertainty when the NNLO corrections to the qq channel are included.
In order to assess the relative size of the contributions from the different colour factors included at NNLO in our calculation, in figure 4 we show the breakdown of the NNLO corrections to the qq channel as functions of p t T and y t . We find that the leading-colour piece, proportional to N 2 c , contributes most significantly, followed by N l N c , which was calculated in [17]. The contributions from all other colour factors are very small.
To further disentangle the light and the heavy quark contributions, in fig 5 we show only the fermionic contributions, omitting the N 2 c colour factor.

Differential distributions for Tevatron
In figure 6 we present differential distributions for top pair production in pp collisions with √ s = 1.96 TeV at LO, NLO and NNLO, along with the corresponding ratios. Given the dominance of the quark-antiquark channel in this scenario, the impact of the NNLO corrections included in our calculation is more important here than in the LHC distributions presented above. In all cases we find good agreement with experimental data, and a significant reduction in the scale uncertainty at NNLO.
As can be seen from the distribution in p t T given in figure (a), in the lower part of the spectrum, NNLO corrections introduce a 10% shift with respect to the NLO prediction. This shift decreases as p t T increases, becoming negligible in the tail of the distribution. In the distribution in |y t | given in figure (b), NNLO corrections to the qq channel introduce an increasing shift of up to 30% with respect to the NLO result. A significant reduction in the scale uncertainty is also observed. In figure (c) we show the inclusive cross section for top pair hadro-production as a function of the invariant mass of the top-antitop system. NNLO corrections, in this case, cause a positive shift over the entire spectrum ranging from 15% near the production threshold to approximately 5% in the tail of the distribution.
It is well known that QCD corrections to the y t distribution in pp collisions cause a shift towards the forward region. This fact can be observed in our result shown in figure (d), where it can also be seen that NNLO corrections modify this shift in comparison to the NLO result. As we will see in the next section, this will have an effect in the NNLO result for the forward backward asymmetry.
In figure 7 we show the breakdown of the NNLO corrections to the y t and m tt distributions into colour factors. The same pattern found in figure 4 can be observed, with the N 2 c colour factor contributing most significantly, followed by N l N c . As before, in figure 8, we omit the leading color factor in order to further examine the size of the different fermionic colour factors.

The forward-backward asymmetry at Tevatron
The forward backward asymmetry in hadronic top quark pair production is an observable that measures the difference between forwardly and backwardly produced top quarks. It is defined as with ∆y tt = y t − yt and ∆y tt ± = θ(±∆y tt ). At leading-order A F B vanishes given the forward-backward symmetry of the Born matrix elements. At higher orders, however, loops and real emissions in the quark-antiquark and quark-gluon channels cause the (anti)top to be more likely produced in the hemisphere defined by the direction of the incoming (anti)quark. Due to this fact, which was first  observed in [53], a non-vanishing A F B in pp collisions is predicted starting at next-toleading order.
Much attention was dedicated to the forward-backward asymmetry after measurements by the CDF collaboration showed a pronounced discrepancy with the Standard Model NLO prediction [5,54]. A recent publication by the D0 collaboration [55] showed measurements that differ from those by the CDF collaboration, and are in agreement with all SM predictions. In [14], the NNLO QCD corrections to A F B were calculated for the first time, yielding a significant reduction in the scale uncertainty and confirming the agreement between the SM prediction and the latest result from D0. Renormalisation and factorisation scales are set equal µ R = µ F = µ and varied as m top /2 ≤ µ ≤ 2m top . Experimental data points from CDF and D0 are taken from the results in the + jets channel presented in [3] and [4] respectively.
Here we present our results for the forward-backward asymmetry at NNLO, computed in the so-called unexpanded form: .

(5.2)
In figure 9 we show A F B as a function of |∆y tt | = |y t − yt| as well as in m tt and p tt T . In the first two cases we find agreement with the D0 measurements, as well as a substantial reduction in the renormalisation and factorisation scale dependence at NNLO. The change from NLO to NNLO in the scale variation bands of the p tt T distribution is rather unusual, with a reduction in the scale dependence in the first bin, and an increase in all others. We compared the results in figure 9 for the differential asymmetries with [14]. Although at NNLO our computation only includes the quark-antiquark channel at NNLO (and it does not include all colour factors), when this channel is dominant, the two results are in agreement. Our results appear to differ at the edges of the distributions. In those regions the top quark pair are produced mostly asymmetrically where we expect the quark-gluon channel, not included in our computation (but included in [14]), to play an important role.

Conclusions
We presented the computation of the virtual-virtual two-parton final state contributions to top pair production in the quark-antiquark channel proportional to the leading colour factor  Figure 9. Differential A F B in: (a) absolute rapidity difference |∆y tt | = |y t − yt|, (b) invariant mass of the tt system m tt , and transverse momentm of the tt system p tt T . The NNLO contributions included are in the quark-antiquark channel only. Renormalisation and factorisation scales are set equal µ R = µ F = µ and varied as m top /2 ≤ µ ≤ 2m top . Experimental data points from CDF and D0 are taken from the results in the + jets channel presented in [54] and [55] respectively. N 2 c . Together with the three and four-parton contributions derived in a previous publication [17], this enabled us to complete the NNLO corrections to top pair hadro-production in the quark-antiquark channel for the phenomenologically most important colour factor.
We derived the subtraction term needed at the two-parton level using the antenna subtraction method extended to deal with massive coloured final states. This required the integration of new massive tree-level four-parton and three-parton one-loop antennae, which we presented here for the first time. We showed that the virtual-virtual counter-term has a similar structure to that found in the massless case [18], and can be expressed in terms of integrated massive dipoles whose pole part is related to the massive infrared singularity operators of [44]. We furthermore demonstrated that the explicit poles of the double virtual matrix elements can be analytically cancelled against those of the corresponding counter-term,

JHEP12(2015)074
providing a crucial check on the correctness of our calculation. In addition, we presented the computation of the heavy quark contributions related to the colour factors N h N c and N h /N c . In particular, we presented our results using analytic expressions for the loop amplitudes in the so-called decoupling scheme, which were not available as such in the literature.
Together with the NNLO corrections proportional to the number of light quark flavours N l obtained in [29], as well as other easily added fermionic colour factors, we implemented the leading-colour and heavy quark NNLO contributions to qq → tt + X completed in this paper, in a fully differential Monte Carlo event generator. This program allowed us to produce several differential distributions for top pair production with NNLO corrections in the quark-antiquark channel for Tevatron and LHC energies. In particular, we obtained distributions in the top quark transverse momentum p t T and rapidity y t , as well as in the invariant mass and rapidity of the top-antitop system m tt and y tt . As expected, we found that at the LHC, the NNLO corrections to the qq initiated process are not phenomenologically significant, except in the very forward and backward regions of the rapidity spectrum. For Tevatron, on the other hand, we found that these corrections are important and drastically reduce the scale dependence in all distributions considered.
In addition, we evaluated the NNLO corrections to the differential forward-backward asymmetries at Tevatron and found good agreement with the results from D0 [55]. In the regions where the contributions from the quark-gluon channel are expected to be small, we also found agreement with [14].

JHEP12(2015)074
with the solid angle given by We can now change variables to the following invariants: The corresponding Jacobian is Note that the energy E 2 is fixed such that 2p 2 · q = E 2 cm + Q 2 and s 12 is not an independent integration variable. It is given by Using the above expressions, the phase space measure in d = 4 − 2 dimensions reads: We now define the following dimensionless variables: and and (B.17) In terms of the dimensionless variables z ij the phase space reads, We use the delta functions to integrate out z 12 and z 34 , and reparametrise the remaining integration variables as where z ± 14 are the roots of∆ 4 viewed as a function of z 14 . The Jacobian for this reparametrisation is dz 13 dz 23 dz 24 = χ 1 (1 − χ 1 )(1 + u 0 − (1 − χ 1 )χ 2 ) (1 + u 0 ) 2 (u 0 + χ 1 ) dχ 1 dχ 2 dχ 3 , (B.21) and the integration region gets mapped to the hypercube 0 ≤ χ i ≤ 1. We finally arrive at . (B.22) In the soft limit x i → 1, after the trivial integration over χ 4 is performed, our phase space parametrisation completely factorises as, Using eqs. (B.23) and (B.24) we can find explicit all order expression for all master integrals that do not involve the invariants s 14 or s 34 . As a boundary condition for the integrals found in the reduction of A 0 q,Qgg we only needed to compute the soft limit of I [4,7,8] , which reads I (sof t) [4,7,8]  In this appendix we collect the seven master integrals found in the IBP reduction of the integrated antenna A 1,lc q,Qg . The Laurent expansion of these integrals can expressed in terms of HPLs with arguments x i or x 0 and GPLs of argument x 0 and weights involving 1/x i .
The phase space associated to the DIS-like kinematics of this antenna, namely p 2 + q → p 1 + p 3 (C.1) with p 2 2 = p 2 3 = 0, p 2 1 = m 2 Q , q 2 = −Q 2 < 0, was derived in [24]. It is given by The inclusive phase space integral, denoted as I [0] , is obtained by performing the y integration in eq. (C.2): With this parametrization of the phase space, the invariants can be expressed as . (C.5) -41 -