NNLO QCD subtraction for top-antitop production in the $q\bar{q}$ channel

We present the computation of the double real and real-virtual contributions to top-antitop pair production in the quark-antiquark channel at leading colour. The $q \bar q \to t \bar{t} g$ amplitudes contributing to the real-virtual part are computed with OpenLoops, and their numerical stability in the soft and collinear regions is found to be sufficiently high to perform a realistic NNLO calculation in double precision. The subtraction terms required at real-real and real-virtual levels are constructed within the antenna subtraction formalism extended to deal with the presence of coloured massive final state particles. We show that those subtraction terms approximate the real-real and real-virtual matrix elements in all their singular limits.


Introduction
Top quark physics has become precision physics at the LHC. Some observables, like the total cross section for tt production, are expected to be measured with accuracies at the percent level. In addition, the ATLAS and CMS collaborations at CERN have reported first measurements of differential observables in top-quark pair production, such as the transverse momentum and rapidity of the tt system [1], its invariant mass [2], and the top quark transverse momentum [3]. Those measurements will allow for a much more detailed probe of the top quark production mechanism than what can obtained from the total cross section. To reliably interpret these data, these precise measurements have to be matched onto equally accurate theoretical predictions. Those can be obtained by computing these hadron collider observables at the next-to-next-to leading order (NNLO) in perturbative QCD. At present, a fully differential NNLO calculation of the cross section for top pair production including all partonic channels is still missing. Intermediate results have recently become available in [4][5][6][7][8][9][10][11][12][13][14][15].
Most notably, the inclusive total hadronic tt production cross section has been presented in [16].
At NNLO, perturbative calculations of collider observables, like jet or heavy quark cross sections and associated kinematical distributions, are typically carried out using partonlevel 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. Towards this ultimate goal for top-pair production observables, in this paper we consider the quark-antiquark initiated channel at leading colour and compute two essential contributions to the NNLO top pair production cross section, namely the double real and real-virtual parts.
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 production process, denoted as double-virtual contributions dσ V V , the one-loop corrections to the (n+1)-parton production process, called real-virtual contribution dσ RV , and the tree-level (n + 2)-parton double real contribution, dσ RR . 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, the real-virtual and virtual-virtual contributions to hadron collider observables involve initial-state collinear singularities that must be absorbed into mass factorisation counter terms. Those are labelled as dσ M F,1 and dσ M F,2 , respectively.
The combination of subprocesses of different particle multiplicity and the consistent cancellation of the respective infrared singularities is one of the major challenges in the construction of NNLO parton-level event generators. In each subprocess, infrared singularities assume a different form: in the virtual corrections they are explicit, while in the real contributions they are implicit and become explicit only after phase space integration. To compute an observable beyond leading order, a regularization procedure is therefore required to extract and cancel the infrared singularities among different partonic channels before those can be implemented in the parton-level event generator. This goal is typically achieved by means of subtraction methods, where all relevant singularities of the matrix elements are subtracted by means of universal auxiliary terms, which are sufficiently simple to be added back after analytic integration over the unresolved phase space. In the past, this approach was successfully applied to various NNLO calculations using sector decomposition [17][18][19][20], q T -subtraction [21], antenna subtraction [22] and most recently with an approach based on sector-improved residue subtraction [12,23].
Two of these methods have been extended to treat massive final state fermions and applied to top pair hadro-production. In [16] the total cross section for inclusive tt production was obtained with the Stripper method [23,24], which combines the FKS subtraction method [25] and sector decomposition [18,19]. Moreover, the antenna subtraction formalism with massive fermions has been applied to the evaluation of the double real contributions to tt production for the pure fermionic processes [4] and for the gluon initiated process gg → ttqq [26]. In this paper, we shall employ the massive extension of antenna subtraction to extract the infrared behaviour of double real and real-virtual NNLO contributions to the qq → tt channel at leading colour.
While the computation of NNLO corrections to observables involving massive par-ticles require the same kind of ingredients as for massless observables, namely real-real, real-virtual and virtual-virtual contributions, the presence of massive fermions in the final state introduces a few simplifications as well as new complications. First, due to the presence of massive final states, the ultraviolet renormalisation procedure of one and two loop amplitudes is more involved than for their massless counterparts. Not only couplings but also mass and wave function ultraviolet renormalisations are required. For all loop amplitudes encountered in this paper, we shall use the ultraviolet regularisation procedure described in [8,10]. Concerning infrared singularities, massive quarks do not give rise to final-state collinear singularities, and the quasi-collinear effects described in [27,28] can be safely ignored for tt production at the LHC. Thus only divergencies associated with soft radiation and with collinear emissions off massless partons require explicit subtraction terms. On the other hand, the non-vanishing parton masses introduce a new scale, which represents a considerable source of complexity both for the final-state kinematics and for the integration of the subtraction terms.
Employing a subtraction method, the NNLO partonic cross section for top-pair production in a given partonic channel (and proportional to a specific colour factor) has the general structure [22] 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 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 antenna subtraction framework employed in this paper, we decompose further the double real subtraction term dσ S NNLO . This term contains distinct pieces corresponding to different limits and different colour-ordered configurations. Some of these pieces ought to be integrated analytically over the unresolved phase space of one particle and combined with the 3-parton final state, while the remaining terms are to be integrated over the unresolved phase space of two particles and combined with the 2parton contributions. This separation amounts to splitting the integrated form of dσ S NNLO as [29][30][31][32] which allows us to rearrange the different terms in eq.(1.1) into the more convenient form In this paper, we shall explicitly construct the antenna subtraction terms dσ S NNLO and dσ T NNLO entering at the four-and three-parton contributions to the NNLO top pair production cross section (1.3) for the quark-antiquark channel at leading colour. The virtual-virtual subtraction term, dσ U NNLO , will be derived elsewhere. Based on the universal factorisation properties of QCD colour-ordered amplitudes, the antenna formalism [4,22,26,27,31,[33][34][35][36][37][38][39][40][41][42] provides a general framework for the construction of subtraction terms that reproduce the singular behaviour of the double real and mixed real-virtual NNLO corrections. Subtraction terms are constructed as products of antenna functions and reduced matrix elements squared with remapped momenta, and the subtraction procedure is based on the colour ordering representation.
The antenna functions capture all the unresolved radiation emitted between a pair of hard partons, referred to as hard radiators. In hadronic collisions, the hard radiators can be initial or final state partons, and therefore three types of antennae must be distinguished: final-final (f-f), initial-final (i-f) and initial-initial (i-i). While NLO subtraction terms only involve three-parton tree-level antennae, at NNLO four-parton tree-level antennae and three-parton antennae are also needed in the double real and real-virtual contributions, respectively. In addition, 3 → 2 and 4 → 2 phase space mappings are required for the reduced matrix elements multiplying the antenna functions in the subtraction terms. Moreover, the analytic integration of the subtraction terms over the appropriate unresolved patch of the phase space requires an exact and Lorentz invariant factorisation of the phase space. Both the mappings and the factorisations are different in f-f, i-f, and i-i configurations. They can all be found for the massive case in [4,27].
The framework outlined above for the construction of NNLO antenna subtraction terms was set up in [31,32,42] in the context of a proof-of-principle implementation of the purely gluonic leading-colour NNLO contributions to di-jet production at hadron colliders. In [31,32], the correctness of the method was checked by showing a complete cancellation of all explicit and implicit infrared divergences that arise in the intermediate steps of the calculation. These results were numerically implemented in the NNLOJET parton-level event generator [43], producing the first NNLO results for hadronic di-jet production. A considerable reduction of the theoretical scale uncertainty was observed, and for the first time double differential distributions in p T and η for inclusive-jet and di-jet NNLO cross sections were presented. Recently, these results have been upgraded to include the full colour dependence in [29].
As outlined above, the goal of this paper is to employ the antenna subtraction method in its extension to the massive case to compute the double real and real-virtual corrections to tt hadro-production in the qq channel. In particular we shall focus on the leading-colour pieces of processes qq → ttgg at tree-level and qq → ttg at one-loop and their corresponding antenna subtraction terms denoted as dσ S NNLO and dσ T NNLO in eq.(1.3) . This will require the use of known phase space factorisations and mappings [4,27,44] and several massive antennae. From those, the three-parton massive antennae are known. New tree-level fourparton and three-parton antennae will be derived here for the first time. The general structure of the subtraction terms remains unchanged with respect to the massless case [22,31,32,42] though, and we shall not repeat it here. We shall instead restrict our presentation to the new elements that are relevant for the present calculation and are related to the presence of massive final states.
Besides antenna subtraction, also the calculation of the 2 → 3 one-loop amplitudes represent a nontrivial ingredient of the 2 → 2 NNLO calculation at hand. Thanks to the recent advent of fully automated NLO tools, such contributions can be in principle computed on a routinely basis. However, the application of NLO tools in the framework of NNLO calculations poses new and still poorly explored challenges. First of all, depending on the employed tool, the numerical character of the new one-loop algorithms might imply a serious CPU speed penalty as compared to analytic approaches. Moreover, the integration of the (subtracted) one-loop contributions over the soft and collinear regions of phase space can lead to serious numerical instabilities. In particular, the well known spurious singularities related to small Gram determinants are inevitably enhanced in the infrared phase space regions, and the resulting loss of numerical accuracy can be strongly enhanced by the large cancellations between matrix elements and subtraction terms. It is thus a priori not clear if automated one-loop generators can guarantee an adequate level of numerical stability and speed for NNLO applications. In this paper we address these issues using the OpenLoops [45] one-loop generator in combination with the Cuttools [46] reduction library, which allows us to study the behaviour of one-loop matrix elements in the deep infrared regime using quadruple precision. As we will show, in spite of the presence of severe instabilities associated with very soft gluon emissions, in the antenna subtraction framework the employed tools turn out to be sufficiently stable to perform a realistic NNLO calculation in double precision. Given the high speed of OpenLoops and the fact that quadruple precision can be avoided almost completely, this guarantees a highly efficient integration of the real-virtual NNLO contributions.
The paper is organised as follows: In section 2, we shall present the cross section for the top-antitop production up to the NLO level. This will enable us to set up the normalisation and present the NLO ingredients required for the computation of the top-pair production cross section at NNLO. The new tree-level four-parton antenna and the three-parton antenna functions required at the double real and real-virtual level of this computation will be presented in sections 3 and 4, respectively. The double real contributions and their subtraction terms are derived in section 5. Sections 6 and 7 contain the real-virtual contributions. Their general structure is presented in section 6 while their computation with OpenLoops is described in section 7. In section 8, we explicitly construct the real-virtual subtraction terms which cancel the explicit infrared poles of the real-virtual contributions and approximate these contributions in all their single unresolved limits. Sections 9 and 10 present various detailed checks on the consistency and numerical stability of the double real and real-virtual subtractions. Finally, section 11 contains our conclusions. In appendix A, the universal single unresolved soft and collinear factors are presented, in appendix B the colour-ordered infrared singularity operators are included, while in appendix C the full expression of the antenna presented in section 3 is given.
2 Top-antitop production in the qq channel at NLO In this section we shall present the main ingredients that are required in the computation of the NLO cross section for tt hadronic production in the qq channel. Besides setting up the notation and the general framework that we will follow throughout the present paper, those NLO contributions will be an essential input for the NNLO mass factorisation counter term dσ MF,1 qq,NNLO which shall be derived explicitly in section 8.

Notation and conventions
To facilitate the reading of our expressions, we shall closely follow the notation in [4,27] for the matrix elements and subtractions terms. In order to identify the colour-ordered subamplitudes with the colour factors that multiply them in the colour decomposition of the full amplitude, we use the following conventions: Different colour strings are separated with double semicolons. A colour string (T a 1 . . . T an ) ij corresponds to . . . ; ; i, a 1 , . . . , a n , j; ; . . . in the argument of the corresponding sub-amplitude. Adjacent partons within one colour string are colour-connected. An antiquark (or an initial state quark) at the end of a colour chain and a like-flavour quark (or initial state antiquark) at the beginning of a different colour chain are also colour-connected, since the two chains merge in the collinear limit where the qq pair clusters into a gluon. When decoupling identities are used, we denote the gluons which are photon-like and only couple to quark lines with the index γ instead of g. In sub-amplitudes where all gluons are photon-like no semicolons are used, since the concept of colour connection in not meaningful. Finally, a hat over the label of a certain parton indicates that it is an initial state particle (for example,1 q is an initial state quark with momentum p 1 ).
Concerning the kinematics, we will use the following definition of invariants: both for massless and massive momenta to make the mass-dependent terms explicitly proportional to m Q . The momenta p i,j,k in eq.(2.1) have to be understood as physical incoming or outgoing momenta with p 0 i,j,k > 0. The above invariants are thus always positive, and crossing transformations have to be accompanied by sign-flips s ij → −s ij , s ijk → −s ijk whenever appropriate.

tt production at LO
The hadronic cross section for tt production at leading order involves two partonic channels, with either a qq pair or a pair of gluons in the initial state. It is given by where H 1 and H 2 are the momenta of the incoming hadrons, p i = ξ i H i , and the sum runs over all quark flavours. Restricting ourselves to the qq initiated process, the leading order partonic cross section takes the form: where, dΦ 2 (p 3 , p 4 ; p 1 , p 2 ) is the 2 → 2 partonic phase space, J 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 the colour-ordered and coupling-stripped tree-level amplitude. It is related to the full amplitude through the (trivial) colour decomposition The normalisation factor is where s is the energy squared in the hadronic center-of-mass frame. Included in this normalisation factor are the flux factor, as well as the sum and average over colour and spin. The constants C( ) andC( ) are defined as: providing the useful relation

tt production at NLO
At the next-to-leading order, three different partonic channels enter: the qq, the gg and the qg channels. The hadronic cross section for tt production at this order is therefore given by where we have used the fact that the partonic cross sections for the qg and theqg are identical due to their invariance under charge conjugation. Restricting ourselves to the qq channel and employing a subtraction method at NLO, the partonic cross section takes the form, (2.9) The three-parton final state contains the real radiation contributions dσ R qq,NLO and their corresponding subtraction term dσ S qq,NLO , whereas the two-parton final state includes the virtual contributions dσ V qq,NLO , the integrated subtraction term 1 dσ S qq,NLO and the NLO mass factorisation counter term dσ MF qq,NLO . The latter is related to the leading order partonic cross section and will be given below.
By grouping the different contributions to the NLO partonic cross section as in eq.(2.9), the difference (dσ R qq,NLO − dσ S qq,NLO ) is numerically well behaved in all regions of the 2 → 3 phase space. It can be integrated numerically in four dimensions. Furthermore, the twoparton contributions (dσ V qq,NLO + 1 dσ S qq,NLO +dσ MF qq,NLO ) are free of poles in the dimensional regulator as we shall demonstrate below.

Real radiation contributions
The real radiation corrections to the qq channel for tt hadronic production are due to the process qq → ttg. The colour decomposition of the corresponding tree-level amplitude is (2.10) Squaring this expression and combining it with the 2 → 3 phase space, the appropriate overall factors and the measurement function, we can write the real radiation contributions as where we have defined in which the gluon is U (1)-like. The normalisation factor N R,qq NLO is given by 13) and the measurement or jet function denoted by J 3 (p 3 , p 4 , p 5 ) guarantees that out of three-parton with momenta p 3 , p 4 and p 5 a final state with a massive heavy quark pair is formed.
The matrix elements squared in eq.(2.11) can become singular when the gluon, whose momentum is denoted by p 5 in the above equation, is either soft or collinear to either of the incoming partons. The antenna subtraction term that reproduces the behaviour of dσ R qq,NLO in those limits is known [27]. It is constructed entirely with products of Atype antennae and reduced matrix elements in final-final, initial-final and initial-initial kinematical configurations (2.14) The reduced matrix elements and measurement functions in the equation above contain redefined momenta that are obtained from the original ones through Lorentz invariant onshell mappings, whose form is different in subtraction terms involving final-final, initial-final and initial-initial antennae. Final state and initial state remapped momenta are denoted with tildes (e.g. p 35 ) and bars (e.g.p 1 ), respectively. In final-final subtraction terms the mappings employed are of the form {p i , p j , p k } → { p ij , p jk } and both redefined momenta are obtained from all three original momenta in the antenna system. Initial-final mappings are of the form {p i , p j , p k } → {p i , p jk }, where p i is an initial state momentum which the mapping rescales, and p jk is obtained from all three momenta in the antenna system. For subtraction terms in initial-initial configurations, the mapping rescales both initial state momenta and performs a Lorentz boost on all remaining final state particles in order to preserve momentum conservation in the reduced matrix elements. The precise definitions of all these mappings can be found, for example, in [27].
The construction of the subleading colour pieces (1/N c ) of eq.(2.14) requires a special procedure, which was explained in [4,27].
The integrated form of the NLO subtraction term dσ S qq,NLO is obtained by factorising the 2 → 3 phase space into an antenna phase space and a reduced 2 → 2 phase space, and integrating the antenna functions A 0 3 in eq.(2.14) inclusively over the antenna phase space. This factorisation, as well as the specific form of the antenna phase space is different for final-final (f-f),initial-final( i-f) and initial-initial (i-i) configurations. It has been derived in the massless case in [22,36] and in the massive case in [27,40]. The integrated forms of the A-type antennae in eq.(2.14) are denoted as A. We shall only make explicit use of their pole parts, which can be entirely written in terms of universal splitting kernels and infrared singularity operators as The colour-ordered infrared singularity operators of the form I (1) ij appearing in the above equation are given in appendix B. The splitting kernel Γ (1) qq (x) in D = 4 − 2 dimensions is given by From these equations, we can express the pole part of the integrated form of eq.(2.14) compactly as Qq ( , s 14 )

Virtual contributions
The virtual contributions denoted by dσ V qq,NLO in eq.(2.9) are due to the process qq → tt at . The colour decomposition of the relevant one loop amplitude is Each of the partial amplitudes can be still decomposed into primitives as (2.21) where N l and N h are respectively the number of light and heavy flavours. Using eqs. (2.20) and (2.21) together with the colour decomposition given in eq.(2.4) for the corresponding tree-level amplitude, we can write the virtual contributions to the tt production cross section in the qq channel as where we have introduced the following compact notation Interestingly, the partial amplitude M 1 4,2 (3 Q ,1 g ,2 g , 4Q) present in eq.(2.20) vanishes when interfered with the tree-level amplitude of eq.(2.4) and it drops out of dσ V qq,NLO . These virtual contributions have been computed in [47][48][49]. Our expressions are in full agreement with those known results.
The matrix elements in eq.(2.22) contain ultraviolet as well as infrared divergences. While the infrared divergences cancel when added to the integrated subtraction terms and mass factorisation terms, the UV poles are removed by renormalisation. For all loop amplitudes throughout this paper we shall follow the renormalisation scheme described in [8,9], in which the heavy quark mass and wave function are renormalised on shell, while the strong coupling constant is renormalised in the MS scheme. In the particular case of the amplitude for the process qq → QQ no mass renormalisation is needed since the corresponding tree-level process does not contain any internal massive propagators. With this simplification, the amplitude is renormalised as where the subscripts (bare) and (ren) stand for bare and renormalised respectively, and the renormalisation constants are given by The explicit infrared pole structure of the UV-renormalised virtual contributions dσ V qq,NLO can be casted in terms of massless and massive colour-ordered infrared singularity operators I (1) ij as, (2.28) As can be seen in the equation above, after UV renormalisation, the remaining infrared poles of the virtual contributions are proportional to the colour factors N c and 1/N c . The absence of infrared poles in the closed-fermion-loop contributions, that is, the contributions proportional to N l and N h , is expected, since the real radiation contributions in eq.(2.11) have no terms proportional to N l or N h . We have cross checked eq.(2.28) against the known universal pole structure of QCD amplitudes with massive fermions [50], and found complete agreement.

The mass factorisation counter term at NLO
The general form of the NLO mass factorisation counter term is related to the leading order cross section through with the kernel Γ (1) ij;kl defined as and Γ (1) ab (z) are Altarelli-Parisi kernels [51]. Applying this expression to the qq initiated channel we find where1 q and2q have momenta x 1 p 1 and x 2 p 2 respectively and Γ (1) qq (x) was given in eq.(2.16).
Combining eqs.(2.18), (2.28) and (2.31), it is straightforward to see that Within the antenna formalism [22], the singular limits of the double real contributions that occur when a pair of colour-connected partons become simultaneously unresolved are captured by tree-level four-parton antenna functions. In general, these four-parton antenna functions are denoted as X 0 4 (i, j, k, l), and depend on the parton momenta p i , p j , p k , p l and the masses of the hard radiators m i and m l in the massive case. They are obtained from ratios of colour-ordered matrix elements squared as where S ijkl,IL denotes a symmetry factor associated with the antenna which accounts both for potential identical particle symmetries and for the presence of more than one antenna in the basic two-parton process. This factor is fixed by demanding that the antennae collapse exactly into the unresolved factors appropriate to each unresolved limit. The flavours of the partons I and L in the two-parton matrix element are determined by the flavour of the two particles that the matrix elements M 0 4 (i, j, k, l) collapses onto when j and k become unresolved. According to the species of partons I and L, antennae can be classified as quark-antiquark, quark-gluon, and gluon-gluon antennae, and depending on whether the hard radiators i and l are in the initial or in the final state, we distinguish between final-final (f-f), initial-final (i-f)and initial-initial (i-i) antennae.
In the context of this paper, one new massive tree-level four-parton antenna is needed to capture the double unresolved behaviour of the real matrix elements squared associated to the partonic channel qq → ttgg in the leading colour component. It is an A-type initialfinal flavour-violating antenna which is denoted as A 0 4 (1 Q , 3 g , 4 g ,2 q ). It is evaluated from the flavour violating tree-level process γ * q → Qgg through the ratio The full expression of this antenna is rather lengthy and it will be left for appendix C. In the remaining part of this section, we shall present the single and double unresolved limits of this antenna A 0 4 (1 Q , 3 g , 4 g ,2 q ). We will start by presenting the double unresolved factors related to its double unresolved limits. The single unresolved factors are well known and are collected in appendix A for completeness. The integrated form of A 0 4 (1 Q , 3 g , 4 g ,2 q ) is presently unknown, it will be obtained by integrating the antenna over the appropriate initial-final antenna phase space, using the techniques developed in [44]. This integrated form will be part of dσ U NNLO given in eq. 1.5 which will be derived elsewhere.

Universal double unresolved factors
When a pair of massless particles becomes simultaneously unresolved, colour-ordered amplitudes squared factorise into a product of a double unresolved factor and a reduced matrix element with two particles less. The form of the double unresolved factor depends crucially on the colour-connection of the unresolved particles: when these are colour-connected a genuine double unresolved factor is obtained, whereas when they are colour-unconnected a product of two single unresolved factors is obtained. In the following we shall present the genuine colour-connected double unresolved factors that we encounter in the unresolved limits of the double real matrix elements squared associated with the partonic process qq → ttgg. These are a massless triple collinear factor corresponding to the triple collinear limit of both final state gluons and one of the initial state fermions, a massive double soft factor, and a massive soft-collinear factor.

Double soft factor of two colour-connected gluons
When two colour-connected gluons j and k become soft between their neighbours i and l an m-particle colour-ordered matrix element factorises as 3) with the double soft current given by [52] Squaring eq.(3.3) and summing over the polarisations of the soft gluons we find with the massive double soft eikonal factor This result converges to the massless colour-ordered double soft factor of [53] in the limit where m i , m l → 0.

Soft-collinear factor in the colour-connected configuration
Soft-collinear singularities occur in those regions of phase space in which a gluon becomes soft and two other massless partons become simultaneously collinear. The factorisation of colour-ordered matrix elements in these limits is different depending on the colour connection of the unresolved particles. When the soft gluon j is colour connected to the collinear particles k and l, the soft-collinear factorisation is given by where P kl→m (z) is one of the single collinear splitting functions in eqs.(A.2-A.4). If either parton k or l are in the initial state, P kl→m (z) will be an initial-final Altarelli-Parisi splitting function.
In the final-final case, the soft-collinear factor Sc i,jkl (m i ) reads whereas in initial-final configurations they are with S ijk (m i , 0) being the single massive soft factor given in eq.(A.10).

Triple collinear factor
In those regions of phase space where three colour-connected massless partons (i, j, l) become collinear, a generic colour-ordered amplitude squared denoted by |M 0 n (. . . , i, j, k, . . .)| 2 factorises as where the three colour-connected final state particles (i, j, k) cluster to form a single parent particle l. The limit is approached in phase space when The triple collinear splitting functions generally depend on the momentum fractions z 1 , z 2 and z 3 , as well as on the invariants s ij , s jk , s ik . The explicit functional form of P ijk→l varies according to the flavours of the three collinear particles as well as on their colour connection. There are two triple collinear splitting functions involving a fermion and two gluons, each of which applies to different colour orderings. In this paper we will need the one corresponding to a colour-ordering of the form . . . ; ; i q , j g , k g , . . ., in which case the splitting function is This triple collinear splitting function corresponds to a configuration in which the three collinear particles are in the final state. However, in the double real corrections for top pair production, only collinear limits of an initial state parton and two final state particles are relevant, given the fact that the tree-level matrix elements contain only two massless final state particles. The initial-final triple collinear splitting functions can always be obtained from their final-final counterparts. For example, the splitting function for the clustering (î, j, k) →l can be related to the final-final case (i, j, k) → l as [54] Pî jk→l (z 1 , z 2 , z 3 , s ij , s ik , s jk ) = where ∆ = 0 if the number of incoming fermions is the same before and after the crossing, and ∆ = 1 otherwise.

Infrared limits of
The four-parton tree-level initial-final massive flavour-violating A-type antenna function denoted by A 0 4 (1 Q , 3 g , 4 g ,2 q ) has the following single and double unresolved limits In this last equation (ang.), stands for angular dependent terms. Those terms arise when a gluon splitting is involved in a collinear limit. In this case, the unresolved single collinear factor is not a spin-averaged Altarelli-Parisi splitting function as given in appendix A but it also involves spin-dependent terms [4].

The massive initial-final antenna
The construction of a subtraction term for the real-virtual corrections to tt production in the qq channel in the leading-colour approximation requires a new initial-final one-loop massive antenna function which we will present in this section.

One-loop antenna functions
Within the antenna formalism, the infrared limits of the real-virtual contributions are captured by three-parton one-loop antennae [22,31]. These are generally denoted as X 1 3 (i, j, k) and they depend on the antenna momenta p i , p j , p k as well as on the masses of the hard radiators in the massive case. In general, these one-loop antenna functions are constructed out of colour-ordered three-parton and two-parton matrix elements as where the tree-level antenna function, denoted by X 0 3 (i, j, k), is given by S ijk,IK denotes the symmetry factor associated with the antenna, which accounts both for potential identical particle symmetries and for the presence of more than one antenna in the basic two-parton process. Initial-final and initial-initial antennae can be obtained from their final-final counterparts by the appropriate crossing of partons to the initial-state. This procedure is straightforward at tree-level but requires some care in the one-loop case, since antennae contain polylogarithms or hypergeometric functions that must be analytically continued to the appropriate kinematical region [37,38].
In any of the three kinematical configurations, the antenna functions can be conveniently decomposed according to their colour factors as follows: 1 In general the sub-antennae have ultraviolet and infrared divergences of explicit and implicit nature. In order to remove the ultraviolet poles, we renormalise the amplitudes in eq.(4.1) following the scheme of [8,9], with the renormalisation constants given in eq.(2.25). We find that the renormalisation prescription of the different sub-antennae is where b 0 = 11/6 and b 0,F = −1/3 are the colour-ordered components of the QCD beta function. We have also defined where M 0 3,1M (i, j, k) is the tree-level amplitude with a mass insertion in the massive propagators. Interestingly, the wave function renormalisation counter terms coming from M 1 3 (i, j, k) cancel against those coming from M 1 2 (I, K), in such a way that the antenna function itself does not require wave function renormalisation.
The antennae that we employ in the real-virtual subtraction terms are renormalised at µ 2 = |s ijk |. To ensure that the matrix elements in the real-virtual contributions and the antennae are renormalised at the same scale, we must substitute After UV renormalisation, one-loop antennae still have explicit and implicit infrared divergences. The structure of the former can be entirely captured by colour-ordered infrared singularity operators; the latter occur when massless partons in the antenna become soft or collinear.

Single unresolved factors at one-loop
The factorisation properties of colour-ordered amplitudes in their soft and collinear limits has been extensively studied in [7,[55][56][57][58][59][60][61][62][63][64][65][66]]. Like at tree-level, the interference of a one-loop amplitude with its tree-level counterpart yields soft eikonal factors and collinear splitting functions in its soft and collinear limits respectively. Those singular factors are also found in the unresolved limits of antennae.
In general, when a gluon becomes soft or a pair of massless partons become collinear, the interference of a one-loop and a tree-level colour-ordered amplitude factorises as where Sing    (4.14) In the following we shall present the explicit form of the singular factors that must be considered in the construction of subtraction terms for the leading-colour real-virtual corrections to top pair production in the qq channel.

Collinear splitting functions
For the partonic process that we are presently considering, i.e. qq → ttg, the splitting function that occurs when the final-state gluon becomes collinear to either of the incoming fermions is P 1 qg←Q (z). In the leading-colour approximation, only the N c part of this splitting function is needed, and it is given by In this equation, z is the momentum fraction carried by the gluon and Pq g→q (z) is the tree-level splitting function whose expression is (4.16)

Massive soft factors
As it occurs at tree-level, when a soft gluon is emitted between massive fermions in the colour chain, the soft factor contains mass dependent terms. While at tree-level the massless soft factor can be obtained from the massive one by setting the massess of the hard radiators to zero, this is no longer the case at the one-loop level: masses are present in the arguments of logarithms that diverge in the massless limit. We must therefore consider separately the soft factors with: (a) two massless hard radiators, (b) one massless and one massive hard radiator, (c) two massive hard radiators. When treating the real-virtual corrections to top pair hadro-production within the leading-colour approximation, only case (b) must be considered. Furthermore, only the N c part of the soft currents and eikonal factors are needed. When a gluon j becomes soft in a primitive amplitude where it is colour-connected to the hard particles i and k, the amplitude factorises as . .), (4.17) where X = lc, l, h, slc, and the tree-level current J µ (p i , p j , p k ) is given by The primitive currents J (1),[X] µ (p i , p j , p k ; m i , m k ) take a different form depending on whether m i and/or m k vanish. These massive soft currents were derived in [7] as tensors in colour space that describe the soft factorisation of full amplitudes rather than of colour-ordered sub-amplitudes. The renormalised colour-ordered currents can be obtained from their results. In the case of one vanishing mass, the leading-colour current that we are presently interested in, reads From the one-loop and tree-level soft currents, massive soft eikonal factors are obtained as S We find with the massive tree-level eikonal factor S ijk (m i , 0) given in appendix A.2.

Infrared properties of A 1,lc
3 (1 Q , 3 g ,2 q ) As mentioned above, in the context of this paper, a new massive antenna is needed to subtract the unresolved infrared limits of the real-virtual contributions related to the partonic process qq → ttg. It is a flavour-violating quark-antiquark antenna denoted by A 1 3 (1 Q , 3 g ,2 q ), which we compute directly in the initial-final kinematics following the definition of eq.(4.1). Working in the leading-colour approximation, only the leading-colour part of the antenna A 1,lc 3 (1 Q , 3 g ,2 q ) needs to be considered. The full expression of this subantenna is too lengthy to be presented in this paper, but its pole part can be compactly written in terms of colour-ordered I (1) ij operators. This pole part will be explicitly needed in section 8 and is given by (4.22) Also the unresolved limits of A 1,lc 3 (1 Q , 3 g ,2 q ) will be required in section 8 in the context of the construction of our real-virtual subtraction terms. They read with the soft and collinear factors defined in eqs.(4.21) and (4.15) respectively. The integrated form of A 1,lc 3 (1 Q , 3 g ,2 q ) is not known at present. It will be part of dσ U NNLO given in eq. 1.5 which will be derived elsewhere.

Double real contributions to qq → tt at leading colour
It is the purpose of this section to present the structure of the double real contributions associated to the tree-level process qq → ttgg at leading colour, and to construct the corresponding subtraction terms. The colour-decomposition of the tree-level amplitude for the partonic process qq → ttgg reads dΦ 4 (p 3 , p 4 , p 5 , p 6 ; p 1 , p 2 ) |M 0 6 (3 Q , i g , j g ,1 q ; ;2q, 4Q)| 2 +|M 0 6 (3 Q , i g ,1 q ; ;2q, j g , 4Q)| 2 + |M 0 6 (3 Q ,1 q ; ;2q, i g , j g , 4Q)| 2 J 2 (p 3 , p 4 , p 5 , p 6 ), where the overall factor 1/2 accounts for the identical gluons in the final state. The normalisation factor is and N qq LO has been given in eq.(2.5). This contribution is singular in several single and double unresolved limits namely The general structure of the double real subtraction terms obtained within the framework of the antenna formalism has been presented in [22,42] in the massless case and extended to the massive case in [4,26]. Without entering into the details of this structure, let us recall that in general, double real antenna subtraction terms, which reproduce the behaviour of the double real contributions in all their single and double unresolved limits, contain five different configurations corresponding to: (a) One unresolved parton (b) Two colour-connected unresolved partons (colour-connected) (c) Two unresolved partons that are not colour-connected but share a common radiator (almost colour-connected) Table 1. Type of contribution to the double real subtraction term dσ S NNLO , together with the integrated form of each term. The unintegrated antenna and soft functions are denoted as X 0 3 , X 0 4 and S while their integrated forms are X 0 3 , X 0 4 and S respectively. M 0 n denotes an n-particle tree-level colour-ordered amplitude.
(d) Two unresolved partons that are well separated from each other in the colour chain (colour-unconnected) (e) Compensation terms for the over subtraction of large angle soft emission.
The antenna content of the subtraction terms for each of these configurations is the same for the final-final, initial-final and initial-initial configurations and it is summarised in Table  1, which is taken from [31]. 2 For the evaluation of the NNLO corrections to heavy quark pair production in the qq channel, the configurations (c) and (e), which always occur together, are not needed and will not be discussed here either. Only (S, a), (S, b) and (S, d) subtraction terms are needed to approximate the double real contributions of eq.(5.2), such that the total subtraction term is given by dΦ 4 (p 3 , p 4 , p 5 , p 6 ; p 1 , p 2 ) 2 (p 3 , p 4 , p ji ) 2 As discussed in [4,29] for example, this content is strictly valid only for leading-colour like double real contributions which involve colour-ordered matrix elements squared. For subleading colour contributions involving interferences of colour-ordered matrix elements, more antenna functions are needed.
All three-parton antennae present in this subtraction term have been derived in unintegrated and in integrated form in [27,36,40]. Furthermore, as can be seen from Table 1, the integrated form of dσ S,a qq,NNLO,N 2 c must be added back at the three-parton level, and it will therefore contribute to the real-virtual counter term dσ T qq,NNLO,N 2 c , which will be presented in section 8.
The (S, b) type subtraction term, denoted by dσ S,b qq,NNLO,N 2 c , takes care of the double unresolved limits of the double real contributions in those subamplitudes in which both final state gluons are colour-connected. It is given by, Two different kinds of structures are involved in this subtraction term: X 0 4 × |M 0 4 | 2 and X 0 3 × X 0 3 × |M 0 4 | 2 . The former subtracts the double unresolved limits while introducing spurious single unresolved singularities, whereas the latter removes these spurious limits ensuring that the four-parton antenna is only active in the double unresolved regions. The four-parton antenna A 0 4 (3 Q , i g , j g ,1 q ) present in this equation appears in a subtraction term for the first time. It was discussed in section 3 together with its infrared limits and its explicit form can be found in appendix C.
As shown in Table 1, the brackets in dσ S,b qq,NNLO,N 2 c should be expanded in order to combine its integrated form with the three and two-parton contributions. The pieces involving products of three-parton antennae, which we denote as dσ S,b 3×3 will require analogous methods as developed in [44] and will be addressed elsewhere.
Finally, the subtraction term of type (S, d), denoted by dσ S,d qq,NNLO,N 2 c , is built out of products of two three-parton antennae and four-parton reduced matrix elements squared. Its role in the partonic process that we are presently considering is to ensure the correct subtraction of the initial-final double collinear limits of the double real contributions given in eq.(5.2). It is given by This subtraction term will be added back to the two-parton counter term dσ U qq,NNLO,N 2 c with both three-parton antennae integrated over their corresponding antenna phase space. We shall not discuss this integration in this paper. In section 9 we will present a series of numerical tests that show that the subtraction term dσ S 6 General structure of the real-virtual contributions to qq → tt at leadingcolour The real-virtual contributions to top-antitop production in the quark-antiquark channel are obtained using the interference of the one-loop and tree-level amplitudes for the partonic process qq → ttg. The colour decomposition of the matrix element reads, where each of the sub-amplitudes has the following decomposition into primitives Interfering the matrix element in eq.(6.1) with the tree-level amplitude in eq.(2.10), combining the result with the phase space and the jet function, and retaining only the terms proportional to N 2 c , we obtain dσ RV qq,NNLO,N 2 × M The leading-colour primitive amplitudes in eq.(6.3) contain ultraviolet poles that must be removed by renormalisation. Following the scheme of [8,9], which was described in section 2.3.2, we renormalise the primitive amplitudes as 2 (p 3 , p 4 , p 5 ). (6.6) As we shall see in section 8, these poles will be canceled by the singly integrated double real subtraction terms and mass factorisation counter terms. The implicit infrared poles, on the other hand, originate from the configurations where the final state gluon becomes soft or collinear to either of the incoming particles. Those will be dealt with the genuine real-virtual subtraction term dσ V S qq,NNLO,N 2 c which will also be constructed in section 8.

Real-virtual contributions to top-antitop production in the quarkantiquark channel with OpenLoops
For the calculation of the matrix elements that enter the real-virtual contributions in eq.(6.3) we employ OpenLoops [45], a fully automated generator of one-loop corrections to Standard Model processes. As discussed in the following, OpenLoops builds Feynman diagrams with a recursive algorithm that allows for a fast and numerically stable evaluation of loop amplitudes. The reduction of amplitudes to scalar integrals can be achieved by interfacing OpenLoops to tensor-integral [67,68] or OPP reduction libraries [46,69,70].
In the context of NNLO calculations, the integration of (subtracted) contributions over soft and collinear regions poses non trivial technical challenges as compared to conventional NLO applications. In particular, the loss of precision resulting from the cancellation between amplitudes and subtraction terms in the soft and collinear regions needs to be compensated by sufficiently high numerical accuracy. However, this is quite challenging since infrared singularities tend to amplify numerical instabilities that arise from spurious singularities (like inverse Gram determinants) in the reduction algorithms. It is thus quite interesting to investigate to which extend automated generators can guarantee an adequate level of numerical stability for NNLO calculations. In this respect OpenLoops has already been shown to be successfully applicable to the calculation of the NNLO corrections to pp → Zγ [71]. In this case, using the q T -subtraction technique [21], it was found that the tensor-reduction library Collier [72], which implements the methods of [67,68,73], is sufficiently stable to perform the entire calculation in double precision. Very recently, OpenLoops was also applied to tt production in association with up to two jets at NLO [74], which is closely related to the present NNLO calculation.
In this work, OpenLoops is used to evaluate the amplitudes for qq → ttg. The interference with the related Born amplitudes, the sums over external colours and helicity states, as well as the ultraviolet renormalisation (6.5) are performed in a fully automated way. The UV-finite but still IR-divergent result is returned in the form of a Laurent series, which must be combined with the corresponding subtraction terms. For consistency with the helicity amplitudes implemented in OpenLoops, the tree matrix elements in eq.(6.6) need to be evaluated in D = 4 dimensions. Tree amplitudes (M 0 ) and loop amplitudes (M 1 ) are expressed as sums over corresponding Feynman diagrams, where the colour factors C (d) associated with individual diagrams are factorised, and the corresponding colour-stripped amplitudes are denoted as A (d) k . All colour structures are reduced to a standard basis {C i }, and the colour information needed to build colour-summed squared matrix elements is encoded in the colour-interference matrix, These colour bookkeeping operations are done only once, using a generic and automated algebraic algorithm, during the generation of the numerical code for a particular process. This approach provides high flexibility in the colour treatment, and the leading-colour approximation used in this paper could be easily implemented via a 1/N c expansion of the colour-interference matrix (7.3). Additionally, in order to obtain the leading colour contribution of the counter-term amplitude, the substitutions C F → N c , C A → 0, T F → 0 are applied to colour factors which are attributed to renormalisation constants. The calculation of colour-stripped loop amplitudes within OpenLoops is based on the representation where the denominators D i = (q + p i ) 2 − m 2 i + iε depend on the loop momentum q, external momenta p i , and internal masses m i . The numerator N (d) (q) corresponds to a particular diagram or to a set of diagrams with the same loop topology. It is expressed as a polynomial of degree R ≤ n in the loop momentum, In contrast to traditional approaches, where the above expressions are constructed via explicit insertion of the Feynman rules, the OpenLoops method consists of a numerical recursion that builds the polynomial coefficients N  4) contributions to the numerator are easily obtained in a process-independent way via so-called R 2 counter terms [75].
For the reduction of amplitudes to scalar integrals, the OpenLoops representation (7.4)-(7.5) allows one to use the tensor-integral or OPP reduction techniques. In the former case, the reduction is performed at the level of process-independent tensor integrals, which are then combined with the corresponding coefficients. In this approach, the Collier library implements systematic expansions in Gram determinants and other kinematic quantities [68], which avoid numerical instabilities due to spurious singularities. In the OPP reduction framework, the reduction is performed at the level of the full integrand in eq.(7.4). This requires multiple evaluations of the numerator function, and using the representation (7.5) in combination with the OpenLoops coefficients, N (d) µ 1 ...µr , renders OPP reduction similarly fast as tensor reduction [45].
In Section 9 we will investigate the numerical stability of the amplitudes in the soft and collinear regions using OpenLoops in combination with the OPP reduction library Cuttools [46]. In this context we will exploit the quadruple precision mode of Cuttools both as a rescue system for matrix elements that are not sufficiently stable in double precision, and for precision tests of the real-virtual cancellations in the deep infrared regime.

Real-virtual subtraction terms
The purpose of the real-virtual counter term dσ T can be safely integrated numerically in four dimensions. The generic antenna content of this counter term has been derived for the massless case in [31,32], and it remains unchanged in the massive case. We will here follow the formalism developed in these references, to which the reader is referred for details.
In general, real-virtual antenna counter terms contain singly integrated double real subtraction terms, NNLO mass factorisation counter terms and genuine real-virtual subtraction terms. For the leading-colour contributions to top pair production in the qq channel the counter term has the following structure and its corresponding antenna subtraction term dσ S kl,NLO . It is given by In the context of this paper, the mass factorisation counter term denoted by dσ MF,1a qq,NNLO,N 2 c is constructed as in eq.(8.4) with dσ R qq,NLO given in eq.(2.11). Retaining only the terms with an overall N 2 c we have The integrated antennae in the equation above have been derived in [27]. Only their pole parts will be needed in the context of this paper. The poles of the flavour violating antenna A 0 q,Qg were given in eq.(2.15), and those of the D-type antennae are given by The pole part of the singly integrated real subtraction term denoted as (S, a) is therefore given by Qq ( , s 24 ) Qq ( , s 13 ) 2 (p 3 , p 4 , p 5 ). (8.12) The measurement function J in eq.(8.12) allows the final state gluon in the reduced five-particle matrix elements squared to become unresolved. The corresponding singular limits of this subtraction term are spurious, since they do not correspond to any physical limits of the real-virtual contribution dσ RV qq,NNLO,N 2 c . Those must therefore be cancelled by other terms in dσ T qq,NNLO,N 2 c . We shall shortly see below that this is indeed the case. Combining eqs.(6.6), (8.12) and (8.6) it can easily seen that eq.(8.8) holds. We shall present these three subtraction contributions separately below. Following the general framework described in [31], in order to subtract the single unresolved limits of the real-virtual contributions given in eq.(6.3) we construct our subtraction terms of the type (VS, a) with one-loop antennae multiplied by reduced tree-level matrix elements and one-loop matrix elements multiplied by tree-level antennae. They read,

Construction of dσ
The three-parton antenna A 1,lc 3 appears here in a subtraction term for the first time. This antenna has been presented together with its singular limits in section 4. Its integration over the antenna phase space will as for dσ S,b,4 require application of the methods presented in [44]. Relabelling the final state momenta, we find Qq ( , s 24 ) + I  The singly integrated subtraction term, 1 dσ S,b 3×3 qq,NNLO,N 2 c , on the other hand, is obtained by integrating the "outer" antennae in eq.(5.6) over the corresponding three-parton antenna phase space. We find 2 (p 35 , p 4 )   In order for eq.(8.14) to be satisfied, dσ VS,b qq,NNLO,N 2 c must be constructed in such a way that its pole part is opposite to the equation above. We must therefore identify the integrated antenna functions that yield the I (1) ij operators and splitting kernels in eq. (8.20). In this case, the integrated antennae that should be employed are initial-final flavour-violating A-type antennae, and it can be seen that eq.(8.14) is satisfied if we write 2 ( p 35 , p 4 ) .
2 (p 3 , p 45 ) . For each singular region we used a series of phase-space samples generated with RAMBO [76] by requiring an increasingly small distance, parametrised in terms of appropriate parameters x k , from the relevant singularity. In the next two sections, we will quantify the level of the real-real and real-virtual cancellations as and respectively. To demonstrate the consistency and stability of the subtractions we will show that the δ RR and δ RV distributions converge to zero in all relevant x k → 0 limits. On the right-hand-side of (9.2) the consistent subtraction of explicit infrared singularities in the numerator and denominator is implicitly understood. Each of the employed samples consists of about 10 4 points with √ŝ = 1 TeV 3 and m Q = 174.3 GeV.

Tests of the double real contributions
We start by discussing infrared cancellations for the double real contribution qq → QQgg in leading colour approximation. To this end we generated 2 → 4 phase space points near all possible single and double unresolved limits. The 2 → 4 tree-level matrix elements in (9.1) have been computed with an in-house Mathematica program based on Qgraf [77] and numerically checked against MadGraph [78] for a few phase space points.

Double soft limits
As shown in Fig.1(a), a double soft phase space point is characterised by the heavy quark pair taking nearly the full energy of the event, and therefore a suitable variable to control the proximity of the events to the singular limit is x = (s − s 34 − 2m 2 Q )/s. In Fig.1 subtraction term is larger than δ RR . The good convergence of the subtraction terms to the double real contributions as the singularity is approached can be seen in the fact that the events accumulate more rapidly near δ RR = 0 as the control variable x is taken to be smaller.

Triple collinear limits
Since we do not subtract collinear limits involving the massive fermions because they are regulated by the large value of m Q , the only types of triple collinear limits that we must consider are of initial-final nature as depicted in Fig.2(a). The control variable in this case is defined as x = s i56 /s, where i = 1, 2. In Fig.2(b) we show how, as we take smaller values of x, i.e. as we get closer in phase space to the singularity of the real radiation matrix element, there is a more rapid accumulation of events around δ RR = 0, signalling again that the approximation is correct. These results correspond to the triple collinear limit p 1 ||p 5 ||p 6 . Similar results are obtained for p 2 ||p 5 ||p 6 .

Soft-collinear limits
As shown in Fig.3(a), soft-collinear limits occur when one of the final state gluons becomes soft and the remaining one becomes collinear to an initial state leg. To probe the softcollinear regions of phase space we generate events with a soft gluon and rotate the final state to make the hard gluon collinear to one of the initial state legs. We employ two control variables x and y. If we consider the limit where gluon (5) is soft and (6) becomes collinear to the incoming leg (1), x is defined as x = (s − s 346 − 2m 2 Q )/s and y is given by y = s 16 /s. As can be seen if Fig.3(b) the convergence of the subtraction term to the partonic double real contribution is once more achieved.

Double collinear limits
Due to the fact that the quasi-collinear limits involving the heavy (anti) quark do not require subtraction, the only double collinear limits in which the double real contributions can diverge are the two simultaneous single collinear limits depicted in Fig.4(a). To control the proximity of the phase space points to the double collinear singularity p 1 ||p 5 , p 2 ||p 6 we employ the variable x = s 15 /s = s 26 /s. As can be seen from Fig.4  results show that behaviour of the double real corrections in their double collinear limits, is correctly described by our subtraction terms. Similar results are obtained for the double collinear limit p 1 ||p 6 , p 2 ||p 5 .

Single soft limits
Single soft limits are characterised by having the three hard final state particles taking nearly the full center-of-mass energy of the event leaving one of the final state gluons with an almost vanishing energy. Consequently, if the soft-gluon momentum is p 5 , we define the control variable as x = (s 346 − s − 2m 2 Q )/s. In Fig.5(b) we show how as the singularity is approached by making x closer to zero, events accumulate more rapidly near δ RR = 0. Analogous results are obtained when the soft-gluon momentum is p 6 .

Final-final single collinear limit
As depicted in Fig.6(a), final-final collinear limits occur when the final state gluons with momentum p 5 and p 6 become collinear. This divergence is approached as the ratio x = s 56 /s gets closer to zero.
As discussed previously in [4,26,39,42], because of the presence of angular correlations between the splitting functions and the reduced matrix elements, in single collinear limits corresponding to the gluon splittings g → gg and g → qq, antenna subtraction terms do not reproduce the behaviour of the real radiation matrix elements in an exact point-bypoint manner but in a two-to-two point manner. This is due to the fact that the angular correlations which spoil the convergence are averaged out when a single collinear phase space point is combined with another single collinear point which differs from the original by a π/2 rotation of the collinear pair around the collinear axis. A thorough discussion of this issue can be found in [42]. In the histogram of Fig.6(b) the aforementioned angular averaging has been performed.

Initial-final single collinear limits
The topology of the single initial-final collinear events is illustrated in Fig.7(a), and the corresponding control variable is defined analogously to the final-final case. There are four different collinear limits in the partonic process qq → QQgg, namely p i ||p j with i and j = 5, 6. Fig. 7(b) contains our results for the limit p 1 ||p 6 , which clearly show that the subtraction terms correctly approximate the double real radiation contributions in this limit. The singularity in Fig. 7(b) is parametrised in terms of x = s 16 /s, and similar histograms are obtained for the other three limits of this kind.

Tests of the real-virtual contributions
In this section we study the cancellation between the real-virtual matrix elements and the corresponding subtraction terms. Due to the lower multiplicity of the 2 → 3 final state and the fact that the heavy quark mass regulates all final-final single collinear limits, the singular structure of the real-virtual contributions is simpler than that of the double real pieces. Indeed, only the soft limit p 5 → 0 and the initial-final collinear limits p i ||p 5 (i = 1, 2) must be considered.
The real-virtual cancellations provide a strong check both of the correctness of the subtraction terms presented in section 8 and of the numerical stability of the OpenLoops amplitudes discussed in section 7. In the vicinity of the soft and collinear singularities matrix elements and subtraction terms are strongly enhanced, and the cancellation can amount to several digits. While this requires augmented numerical accuracy in the unsubtracted amplitudes, numerical instabilities related to Gram determinants can be strongly amplified in the vicinity of the singularities. It is thus crucial to prevent that the infrared cancellations are spoiled by numerical instabilities of the amplitudes. To this end, Open-Loops implements an instability trigger, which monitors the numerical accuracy of the results by means of a scaling test. The amplitudes are evaluated a second time by rescaling all dimensionful input parameters by a constant factor ξ, and the output is rescaled back by a factor ξ −d depending on its mass dimension d. The agreement with the original 3 t matrix element serves as an accuracy estimate, and phase-space points that are not sufficiently stable are automatically reevaluated with a rescue system. Results presented in the following have been obtained with Cuttools as a reduction back end of OpenLoops, using the quadruple-precision mode of Cuttools as a rescue system for unstable points. Matrix elements are first evaluated in double precision and are reevaluated in quadruple precision if their estimated double-precision accuracy is less than 3 correct digits or smaller than the observed cancellation δ RV with the subtraction term. The stability of the quadruple-precision output is assessed with an additional scaling test. Due to the fact that the scaling test tends to overestimate the accuracy, following a universal distribution, one must demand for an accuracy which is higher than the cancellation by about the width of this distribution. For calibration we determine the width from double precision scalings, using a quadruple precision result as reference point, finding a width of around one decimal digit. If needed, the accuracy estimate can be improved using multiple scalings. For what concerns the numerical stability of the matrix elements, in the collinear region it turns out that for the depicted values of the control variable, double precision provides sufficient stability (in the sense of the criterion described above) for the vast majority of the phase space points. This also holds in the soft regions with x = 10 −3 and x = 10 −4 . Starting at the soft sample with x = 10 −5 , a sizable fraction of the matrix elements must be evaluated in quadruple precision. However, it should be pointed out that this deep infrared region (x = 10 −5 corresponds to a gluon energy around 5 MeV) is not relevant for physical applications based on antenna subtraction. In fact, as will be shown in section 10, double precision results are sufficiently stable to obtain integrated cross sections with permil level accuracy.
Detailed findings on the numerical stability and the reliability of the trigger system are summarised in Table 2. The trigger system to detect unstable points from scalings can lead to false positive results, meaning that points will be evaluated in quadruple precision although they were actually stable enough. This is a side effect of avoiding false negative results, meaning points which are regarded as stable although they are not. Note that in the x = 10 −6 soft sample even quadruple precision is no more enough to observe full cancellation for all points, and O(5%) of the points are tagged as unstable. This shows in the tail of the corresponding distribution in Figure 8(b), where the two bins around x = 10 −6.5 are populated only by unstable points.
10 Stability of the integration over the three-particle phase space As a further and more realistic test of the stability of the real-virtual matrix elements and of the related subtraction terms we have integrated the difference dσ RV qq,NNLO,N 2 c −dσ T qq,NNLO,N 2 c inclusively over the three-particle phase space employing a parton level event generator. In this integration, we impose a technical cut on the gluon p T using the control variable y cut = p g T / √ŝ , in such a way that no events are generated too close to the soft and collinear singularities. Naturally, since the entire phase space ought to be covered in the integration, y cut must be taken small. While the unsubtracted dσ RV contribution would lead to a logarithmic divergence in the limit y cut → 0, the subtraction term guarantees  Table 2. For the samples of phase space points of Figs. 8(b) and 9(b) the fraction of points is shown which are are unstable in double precision ("unstable"), meaning that the accuracy is not high enough to observe full cancellation between matrix element and subtraction term. "triggered" is the fraction of points which is detected as unstable by the trigger system described in the text, and subsequently evaluated in quadruple precision, and "false negative" is the fraction of points which are unstable, but not triggered. a smooth convergence at small y cut . In practice the integral should reach a plateau for a sufficiently small value of the cut, y max cut , i.e. for any y cut < y max cut the integral of dσ RV − dσ T should remain stable within Monte Carlo integration errors. This is clearly confirmed in for pp → tt as a function of y cut . For both the NNLO real-virtual subtracted contributions and the LO normalisation we used √ s = 7 TeV, m t = 174.3 GeV and set the renormalisation and factorisation scales to µ R = µ F = m t . We employed the MSTW2008nnlo90cl and MSTW2008lo90cl PDF sets for the NNLO and the LO contributions respectively. The high stability of the integration results for values of y cut below y max cut ∼ 10 −3 provides solid evidence of the correctness of the real-virtual subtraction terms of eq.(8.23). Moreover, using OpenLoops in combination with Cuttools, it turns out that the stability plateau is reached before encountering significant instabilities in double precision. For y cut = 10 −3 (10 −4 ) we find that only 1 out of 10 5 (10 4 ) events requires a quadruple precision reevaluation. This allows for a highly efficient evaluation of the real-virtual contributions based on double precision for the vast majority of the phase space points.

Summary and outlook
In this paper, we presented the double real and real-virtual NNLO contributions to hadronic tt production in the quark-antiquark annihilation channel. The computation is performed in leading colour approximation using the antenna subtraction method, which was extended to deal with the presence of a massive fermion pair in the final state. The real-real subtraction terms, presented in section 5, approximate the corresponding 2 → 4 tree matrix elements in all single and double unresolved limits, while the real-virtual subtraction terms, presented in section 8, remove the explicit infrared poles present in the 2 → 3 one-loop matrix elements, as well as the implicit singularities that occur in the soft and collinear limits. The relevant new tree-level four-parton and three-parton massive initial-final antennae functions, together with their unresolved counterparts, have been derived in sections 3 and 4.
The correctness of the subtraction and its numerical stability have been demonstrated with detailed cancellation checks in section 9. To this end, the convergence of the subtracted real-real and real-virtual contributions was studied by means of event samples generated in several phase space slices with increasingly small distance from all relevant single and double-unresolved limits.
To compute the one-loop qq → ttg real-virtual contributions we used OpenLoops in combination with the Cuttools implementation of OPP reduction. This provides interesting insights into the potential benefits of new automated one-loop generators in the framework of NNLO calculations. While the high CPU speed of OpenLoops represents an obvious attractive feature, numerical instabilities could represent a very serious issue for NNLO applications. In fact, while the strong cancellations between one-loop amplitudes and related subtraction terms call for augmented numerical accuracy in the soft and collinear regions, the typical Gram-determinant instabilities of one-loop amplitudes tend to be strongly enhanced in the infrared regions. It is thus important to make sure that the infrared subtractions are not spoiled by numerical instabilities of the one-loop matrix elements. To this end, using scaling tests as well as the quadruple precision mode of the Cuttools library, we performed detailed studies of the interplay between one-loop instabilities and infrared cancellations. On the one hand, it turns out that quadruple precision is essential (and at some point even insufficient) to avoid excessive numerical instabilities in the deep infrared regime. On the other hand, we found that such instabilities arise only at very small gluon energies and are essentially irrelevant for an NNLO calculation based on antenna subtraction. In particular, using a realistic infrared cut-off, one-loop amplitudes in double precision turn out to be sufficiently stable for the vast majority (more than 99.99%) of the phase space points. The fact that quadruple precision can be avoided almost completely implies a drastic efficiency improvement for the integration of the real-virtual NNLO contributions.
In order to complete the NNLO corrections to top-antitop production in the quarkantiquark channel at leading colour, the 2-parton contributions dσ VV NNLO and its corresponding counterterm dσ U NNLO need to be added to the 3 and 4-parton contributions ( dσ RV NNLO and dσ RR NNLO ) and their corresponding subtraction term dσ T NNLO and dσ S NNLO derived in this paper. The 2-loop contributions participating in dσ VV NNLO are known. However dσ U NNLO is presently unknown. In it, essential unknown ingredients are the integrated 4-parton tree level antenna A 0 4 (1 Q , 3 g , 4 g ,2 q ) and the integrated one-loop antenna A 1,lc 3 (1 Q , 3 g ,2 q ) which have been presented in unintegrated form together with their infrared limits for the first time in this paper.
The results presented in this paper constitute a major step towards a complete NNLO calculation, based on antenna subtraction, of top-pair production in the quark-antiquark channel. Our final goal is then the construction of an NNLO parton-level event generator for the two, three and four partonic contributions, which will be applicable to any fully differential observable at hadron colliders.

A.1 The collinear splitting functions
In this paper we have considered the collinear radiation emitted from a massive fermion to be regulated by the mass of this fermion. Consequently, we shall here restrict ourselves to present the usual massless Altarelli-Parisi splitting functions arising in collinear configurations involving only massless partons.
When a pair of massless final state particles i and j with momentum p i and p j become collinear and cluster into a parent parton of momentum p k = p i + p j the kinematics of the collinear configuration can be described as where z is the momentum fraction of one of the partons in the collinear pair. The specific form of the splitting function depends on the species of partons i and j. There are three different splitting functions, corresponding to the three possible final-final parton-parton splittings. In conventional dimensional regularisation, they are given by When one of the collinear particles is in the initial state, the kinematics of the collinear limit can be described as 5) and the four splitting functions corresponding to the four different parton-parton splittings read The additional factors (1− ) and 1/(1− ) account for the different number of polarizations of quark and gluons in the cases in which the particle entering the hard processes changes its type. The antiquark-gluon splitting functions are identical to the quark-gluon ones due to the invariance of the splitting under charge conjugation. In this paper, only the latter splitting functions arising in initial-final collinear configurations are employed.

A.2 The massive soft eikonal factor
When a gluon with momentum p j becomes soft in a colour-ordered tree-level amplitude where it is colour connected to partons i and k with masses m i and m k respectively, the associated soft factor is given by [27,40] S ijk (m i , m k ) = 2s ik s ij s jk − 2m When m i = m k = 0 this factor reduces to the usual massless soft eikonal factor.

B Appendix B: Colour-ordered infrared singularity operators
The explicit pole structure of colour-ordered matrix elements can be written in terms of colour-ordered infrared singularity operators I (1) ij . Within the antenna subtraction method, the pole part of antennae as well as that of integrated tree-level three-parton antennae can be also captured by these operators.
If only massless particles are involved, the following set of operators is sufficient (in addition to the splitting kernels Γ (1) ij (x)) to express the pole structure of a QCD amplitude as well as that of a one-particle inclusive integral of a tree-level amplitude [22,36]: