Top-quark pair production at next-to-next-to-leading order QCD in electron positron collisions

We set up a formalism, within the antenna subtraction framework, for computing the production of a massive quark-antiquark pair in electron positron collisions at next-to-next-to-leading order in the coupling αs of quantum chromodynamics at the differential level. Our formalism applies to the calculation of any infrared-safe observable. We apply this set-up to the production of top-quark top antiquark pairs in the continuum. We compute the production cross section and several distributions. We determine, in particular, the top-quark forward-backward asymmetry at order αs2. Our result agrees with previous computations of this observable.


Introduction
The exploration of the production of top-quark top-antiquark (tt) pairs and their decays is among the core physics issues at future linear or circular electron-positron colliders [1][2][3].
Simulation studies indicate that measurements of the reaction e − e + → tt in the threshold region and in the continuum allow to precisely determine a number of key observables associated with the top quark, including its mass, its width, its Yukawa coupling to the 125 GeV Higgs resonance, and its electroweak neutral current couplings (cf., for instance, [4,5] and references therein). Needless to say, precise predictions are required, too, on the theoretical side. A large effort has been made to investigate tt production at threshold. At present the threshold cross section is known at next-to-next-to-next-to-leading order QCD [6]. As far as the production of tt, or more general, the production of a heavy quark-antiquark pair (QQ) in the continuum is concerned, differential predictions at next-to-leading order (NLO) QCD have been known for a long time for QQ [7] and QQ + jet [8][9][10][11][12][13] final states. Also the NLO electroweak corrections are known [14][15][16][17]. Off-shell tt production and decay including non-resonant and interference contributions at NLO QCD was investigated in [18]. The total QQ cross section σ QQ was computed to order α 2 s (NNLO) and order α 3 s in [19][20][21][22] JHEP12(2016)098 and [23], respectively, using approximations as far as the dependence of σ QQ on the mass of Q is concerned. (A calculation of e − e + → γ * → QQ with full quark-mass dependence of σ QQ was made in [24].) A computation of the cross section and of differential distributions for tt production at order α 2 s with full top-mass dependence was reported in [25,26]. In this paper we set up a formalism for calculating the electroweak production of a massive quark-antiquark pair, e − (p 1 )e + (p 2 ) → γ * , Z * (q) → Q(k 1 )Q(k 2 ) + X , (1.1) at order α 2 s and to lowest order in the electroweak couplings within the antenna subtraction framework and apply it to the production of top-quark pairs. Our approach is fully differential and applies to any infrared-finite observable.
Antenna subtraction is a method for handling infrared (IR) divergences, that is, soft and collinear divergences in higher order QCD calculations [27][28][29][30]. The general features of the method at NNLO QCD were developed in ref. [29]. For QCD processes with massive quarks the antenna subtraction terms at NLO were determined in refs. [31,32]. As to applications of this method to hadronic tt production, partial results were obtained in refs. [33][34][35][36][37]. For the computation of the reaction (1.1) at the differential level we use the unintegrated and integrated NNLO real radiation antenna subtraction terms and the NNLO real-virtual antenna functions worked out in [38,39] and [24], respectively. We recall that alternative methods for handling IR divergences have been successfully applied to NNLO QCD processes involving top quarks. The method of [40,41] was used in the computation of the hadronic tt production cross section [42] and of differential distributions [43]. The results obtained by [25,26] for (1.1) are based on a NNLO generalization of a phase-space slicing method [44,45].
This paper is organized as follows. We recapitulate in the next section the calculation of the differential cross section of (1.1) at order α s using the antenna subtraction method. Section 3 contains a detailed exposition of how to compute within this framework the differential cross section and distributions of IR-safe observables at order α 2 s . In section 4 we apply this formalism to top-quark pair production above the tt threshold. We compute the total tt cross section, a number of differential distributions, and the top-quark forwardbackward asymmetry at order α 2 s . We compare also with existing results. We conclude in section 5. Appendix A contains details about the momentum mappings in the three-and four-particle phase spaces with massive quarks that are required for the antenna subtraction terms at NLO and NNLO QCD.

The differential cross section at LO and NLO QCD
For completeness and for setting up our notation we outline in this section the computation of the differential cross section for e − e + → QQX at order α s within the antenna subtraction method. Here Q denotes a massive quark, for instance, the b or t quark. We work in QCD with n f massless quarks q and one massive quark Q. All matrix elements in this and in the following section 3 refer to renormalized matrix elements. We define the mass of Q, denoted by m Q , in the on-shell scheme while the QCD coupling α s is defined in the JHEP12(2016)098 MS scheme. Dimensional regularization is used to handle IR singularities that appear in intermediate steps of our calculation.

LO QCD
To zeroth order in α s we consider e − (p 1 ) e + (p 2 ) → γ * , Z * (q) → Q(k 1 )Q(k 2 ) , (2.1) where Q denotes a massive quark. The corresponding leading-order (LO) (differential) cross section for unpolarized e − e + collisions is given by where the color-stripped two-parton Born amplitude M 0 2 (1 Q , 2Q) is defined in eq. (2.10) below. We use here and below, as ref. [29], symbolic labels i X in order to display the type X and the four-momentum k i of a final-state parton in the matrix elements. For instance, 1 Q denotes a massive quark with momentum k 1 . Here and in the following, summation over the spins and colors of the partons in the final state is implicit. The factor 1/8s is the product of the spin averaging factor for the initial state and the flux factor, the variable s = (p 1 + p 2 ) 2 , and N c = 3. In general the jet function or measurement function, which must be infrared-safe, is denoted by F (m) n (k j ). It refers to a n-jet observable constructed out of a pair of Q,Q and (m − 2) massless partons in the final state. Here and in the following section we put for definiteness n = 2, but we emphasize that our analysis applies to any infrared-safe observable. The m-particle phase space measure dΦ m in D space-time dimensions is where µ is a mass scale.

NLO QCD
The computation of the NLO QCD correction dσ 1 to (2.2) involves the interference of the Born and one-loop amplitude of (2.1) and the squared Born amplitude of the three-parton final state We recall that in a subtraction scheme for handling the IR divergences, the NLO correction to the LO cross section or to a differential distribution can be written as follows: where = (4 − D)/2 is the parameter of dimensional regularization and the subscripts Φ n denote n-particle phase-space integrals. The second term in the first and second square JHEP12(2016)098 bracket of (2.5) is the unintegrated and integrated subtraction term that renders the difference, respectively sum of the terms in the square brackets finite in D = 4 dimensions. Throughout this paper, the symbol n indicates the analytic integration over the phase space of n unresolved partons in D = 4 dimensions. Within the antenna framework, the NLO subtraction terms required in (2.5) were computed in [31]. The NLO real and virtual corrections to the LO differential cross section, dσ R QQg and dσ V QQ , are given by and γ E = 0.57721 . . . denotes the Euler-Mascheroni constant. As before summation over the spins and colors is understood. We use the following shorthand notation for the interference of the tree-level and one-loop two-parton amplitude: In the formulas (2.2), (2.6), and (2.9) we have introduced color stripped partial amplitudes M 0 2 , M 1 2 , and M 0 3 where QCD coupling factors are taken out, but electroweak couplings are included. These quantities are related to the tree-level and the renormalized one-loop matrix elements of e − e + → QQ and the tree-level matrix element of e − e + → QQg. For reference in the next section we give here the expansion to NNLO QCD of the matrix elements of these processes: where i 1 (i 2 ) denotes the color index of the heavy quark (antiquark), a 3 is the color index of the gluon, g s = √ 4πα s , C F = (N 2 c − 1)/(2N c ), and the generators of SU(3) c are normalized according to tr(T a T b ) = T R δ ab with T R = 1/2. The number of massless quarks is denoted JHEP12(2016)098 by n f . The renormalized two-loop two-parton amplitude M 2 2 , which can be decomposed into different color structures, and the renormalized 1-loop three-parton amplitude in the square bracket of (2.11) are required in the next section. The labels 'lc' and 'sc' in (2.11) refer to leading and subleading color, respectively. The terms M 1,f 3 and M 1,F 3 are the contributions from the massless and massive quark loop, respectively, that enter via the wave-function renormalization of the external gluon. The term M 1,tr 3 denotes the quark triangle contributions, where the axial current couples to quark triangles, summed over all quark flavors (u i , d i ), which disintegrate into a real and virtual gluon that splits into QQ. This term, which is ultraviolet-and infrared-finite, involves weak couplings of q = Q and constitutes in this sense a non-universal correction to the leading-order QQ cross section.
Let's proceed with the discussion of the NLO cross section. The squared Born matrix element M 0 3 2 of the real radiation correction (2.6) diverges when the gluon momentum k 3 becomes soft. Within the antenna method this singularity is regularized by constructing a subtraction term that coincides with (2.6) in this singular limit. The subtraction term and its integrated form (integrated over the phase space of the unresolved gluon) are: where µ is the mass parameter of dimensional regularization, and 14) The three-parton tree-level massive quark-antiquark antenna function A 0 3 and its integrated counterpart A 0 3 were derived in [31,32]. The integrated antenna function A 0 3 contains an explicit IR pole ∝ 1/ that cancels the corresponding IR pole in dσ V QQ . In (2.12) the matrix element M 0 2 and the measurement function F are evaluated with redefined onshell momenta k 13 , k 32 that are obtained from k 1 , k 2 , k 3 by an appropriate phase-space mapping [33]. A method to construct k 13 , k 32 is given in appendix A.1.

The differential cross section at NNLO QCD
The second-order term dσ 2 in the expansion in powers of α s of the differential cross section of (1.1), dσ = dσ LO + dσ 1 + dσ 2 + O(α 3 s ), receives the following contributions: i) the double virtual correction dσ V V NNLO associated with the matrix element of e − e + → QQ to order α 2 s (i.e., 2-loop times Born and 1-loop squared), ii) the real-virtual cross section dσ RV NNLO associated with the matrix element of e − e + → QQg to order α 2 s (1-loop times Born), iii) the double real contribution dσ RR NNLO associated with the squared Born amplitudes e − e + → QQgg, e − e + → QQqq (where q denotes a massless quark). Above the 4Q threshold,

JHEP12(2016)098
e − e + → QQQQ contributes, too. The latter contribution is IR finite and is of no concern for the purpose of this section.
Apart from the QQQQ contribution, the terms i), ii), iii) are IR divergent. Within the subtraction method the second order correction dσ 2 , where the different pieces are separately finite, is constructed schematically as follows: The integrands dσ S NNLO and dσ T NNLO denote the double-real subtraction terms (for QQqq and QQgg) and the real-virtual subtraction term, respectively. We discuss in turn the various terms in (3.1) in some detail.

Double real-radiation corrections
In this subsection we discuss how to compute the first term on the right-hand side of the first line in (3.1) with the antenna subtraction method.
The QQqq final state. First we consider the reaction where q denotes a massless quark. The corresponding tree-level amplitude, decomposed into color-stripped subamplitudes with the QCD coupling factored out, is given by The color indices of the quarks and antiquarks are labeled by i 1 , . . . , i 4 . The matrix element M Q 4 (M q 4 ) corresponds to the subamplitude where the massless (massive) quark-antiquark pair is produced by the splitting of the virtual gluon radiated off one of the quarks produced by the virtual photon or Z boson. The matrix element (3.3) yields the unsubtracted differential cross section summed over colors and summed/averaged over all spins: where the sum is over all n f massless quark flavors and The second and third term in the second line of eq. (3.4) contain the electroweak couplings of the massless quarks q. Thus these terms are non-universal QCD corrections to the leading-order differential QQ cross section. Moreover, these two terms do not become

JHEP12(2016)098
infrared singular in the four-parton phase space. Only the first term in the second line of eq. (3.4) requires subtraction. We define a subtracted differential cross section by subtracting three terms from (3.4).
By construction the sum of the three subtraction terms is such that it coincides with (3.4) in all single and double unresolved limits, i.e., when the massless quarks become collinear and/or soft. Thus, (3.6) is free of IR divergences and can be integrated over the four-parton phase space numerically in D = 4 dimensions. The term dσ S,a,QQqq NNLO subtracts the singularities associated with the single unresolved configurations from the first term in the curly bracket of (3.4). Within the antenna method, it is given by The quark-gluon antenna E 0 3 with a massive radiator (anti)quark is given in [31]. The colorstripped tree-level QQg matrix element M 0 3 is defined in (2.11). The momenta k ij and k jk are redefined on-shell momenta, constructed from linear combinations of the momenta k i , k j and k k [33,46]. The 3 → 2 mappings must be such that they define remapped on-shell momenta that have the correct soft and collinear limits. In appendix A we discuss the 3 → 2 and 4 → 2 mappings that are required for the double real-radiation subtraction terms. A numerical method to construct the mapped momenta that appear in (3.7) is given in appendix A.1.
The subtraction term for removing the singularities of (3.4) due to the double unresolved configuration, where both q andq become soft, is The antenna function B 0 4 is given in [38] and the tree-level matrix element M 0 2 is defined in (2.10). The momenta k ikl and k jkl are linear combinations of the momenta k i , k j , k k , k l obtained from a 4 → 2 mapping, cf. appendix A.2.
The function B 0 4 develops singularities in the single unresolved limits that are subtracted by the additional term

JHEP12(2016)098
The arguments of the antenna functions A 0 3 are mapped momenta obtained by the 3 → 2 mappings described in appendix A.1. The arguments of the Born matrix elements and of the measurement functions in (3.9) are obtained by two consecutive 3 → 2 mappings described in eqs. (A.16) and (A.17) of appendix A.2. These two iterated 3 → 2 mappings are necessary for being able to perform the integration over the antenna phase space of the unresolved parton in analytic fashion and obtain the integrated subtraction term defined in (3.23) below.
Eq. (3.6) is not yet the appropriate expression for numerical evaluation. In the antenna framework there is a subtlety associated with angular correlations [33,[47][48][49]. The gluon radiated off a Q orQ that splits into qq leads to angular correlations in the unsubtracted squared matrix element. However, the type-a subtraction term (3.7) that was constructed to take care of the single unresolved limit of the squared matrix element when the q andq become collinear, is composed of products of spin-averaged three parton antenna functions and three-parton Born matrix elements. That is, the type-a subtraction term does not contain these angular terms and, therefore, does not have the same local singular behaviour as the unsubtracted squared matrix element. The four-parton antenna function B 0 4 and thus the subtraction term (3.8), which takes care of the double unresolved limit, contains these angular correlation terms, but the subtraction term (3.9), which ensures that the complete subtraction term has no singularities in the single unresolved region, does not. However, these angular correlations in the unsubtracted squared matrix element and in the subtraction term (3.8) are averaged out after integration over the azimuthal angle φ between the spatial parts of the light-like momenta k 3 , k 4 and the collinear direction k = k 3 + k 4 . It can be shown [47] that the functional dependence on φ of the squared matrix element in the collinear limit is proportional to cos(2φ + α). This suggests [33,[47][48][49] that the angular correlations can be averaged out by combining, for each final-state momentum configuration, two points in phase space with azimuthal angles φ and φ + π/2. Thus we evaluate M Q 4 2 in (3.4) and the subtraction term (3.8) for each set of momenta k 1 , k 2 , k 3 , k 4 also for k 1 , k 2 , k 3r , k 4r and take the average. The 4-momenta k 3r , k 4r are obtained by rotating the spatial parts of k 3 , k 4 by an angle π/2 around the collinear axis k = k 3 + k 4 . By sampling the phase space in regions where k 3 · k 4 /s 10 −8 we checked that this procedure provides a subtraction term that is a very good approximation to the squared matrix element in the single and double unresolved limits.
The QQgg final state. Next, we consider the reaction (3.10) The corresponding tree-level matrix element can be decomposed into color-ordered subamplitudes as follows:

JHEP12(2016)098
The unsubtracted differential cross section, summed over all colors and summed/averaged over all spins, is given by where the subleading color term is The factor 1/2 in (3.12) is due to Bose symmetry and the factor N is defined in (3.5).
In the subleading color term M sc both gluons are photon-like, i.e., no non-abelian gluon vertices are involved. Hence, when the two gluons become collinear, this term does not become singular.
In analogy to (3.6) we define a subtracted differential cross section by subtracting three terms from (3.12): (3.14) It is by construction free of IR divergences and can be integrated over the four-parton phase space numerically in D = 4 dimensions. As in the QQqq case dσ S, a QQgg NNLO and dσ S, b,2 QQgg NNLO cover the singularities of (3.12) due to single-unresolved and double-unresolved configurations, respectively, that is, when the gluons become collinear and/or soft. The term dσ S,b,1,QQgg NNLO subtracts the singularities of dσ S, b,2 QQgg NNLO in the single unresolved limits.
Within the antenna method, these three subtraction terms are given by dσ S,a,QQgg The tree-level massive quark gluon antenna function d 0 3 is given in [31]. The four-parton QQgg antenna functions A 0 4 andÃ 0 4 , which were derived in [39], govern the color-ordered and non-ordered (photon-like) emission between a pair of massive radiator quarks, respectively. The mapped momenta denoted by a tilde and double tilde in (3.15)-(3.17) are obtained from 3 → 2, 4 → 2, and two iterated 3 → 2 mappings, respectively, in completely analogous fashion as in the QQqq case; cf. appendix A.
The remarks on the angular correlations due to gluon splitting made in the second paragraph below eq. (3.9) apply also here, where the angular correlations are due to g → gg.
These correlations are present only in the leading-color part of the unsubtracted differential cross section (3.12) and in the leading-color part of the subtraction term (3.16). In analogy to the QQqq case we evaluate these leading-color terms for each set of momenta k 1 , k 2 , k 3 , k 4 also for k 1 , k 2 , k 3r , k 4r and take the average. We sampled the phase space for this final state in regions where k 3 · k 4 /s 10 −8 and checked that the resulting subtraction term is a very good approximation to the squared matrix element in all single and double unresolved limits.

Real-virtual corrections
In this subsection we outline how to compute the order α 2 s contribution of the QQg final state to the differential massive quark-pair production cross section with the antenna subtraction method; that is, the second term on the right-hand side of the first line of (3.1).

JHEP12(2016)098
Unsubtracted real-virtual cross section. This contribution involves the interference of the tree-level and one-loop QQg final-state amplitude. Using the conventions of (2.11) the unsubtracted O(α 2 s ) correction to the cross section, summed over colors and summed/averaged over spins, is given by The factors C( ) and N are defined in (2.8) and (3.5), respectively, and we have introduced the shorthand notation with X ∈ {lc, sc, f, F, tr}. We recall that M 1 3 is the renormalized one-loop amplitude. The analytic computation of (3.18), which was first performed in [8][9][10][11][12][13], is standard by now. We recall from section 2.2 that the triangle term δM 1,tr 3 , which was analyzed first in [50], is an IR finite and non-universal QCD correction. The other contributions to the unsubtracted cross section contain explicit IR poles (single and double poles in 1/ ). In addition, the phase-space integration of (3.18) in the region where the external gluon becomes soft, leads to additional IR singularities. Both types of singularities must be subtracted with appropriate terms in order that the integration over the three-parton phase space can be performed numerically in four dimensions.
Subtraction of explicit infrared poles. The explicit IR poles in (3.18) are removed by adding the subtraction terms (3.7) and (3.15), integrated over the phase-space of one unresolved parton: The antenna functions D 0 3 and E 0 3 , which are the integrated versions of the tree-level antenna functions d 0 3 and E 0 3 , respectively, are given in [31,32]. The poles in of these functions and of A 0 3 cancel the explicit IR poles in (3.18). The kinematic invariants that appear in the arguments of these functions are

JHEP12(2016)098
One-loop single-unresolved subtraction term. The singular behavior of (3.18) is mimicked in the limit where the external gluon becomes unresolved by the following subtraction term: The massive one-loop antenna functions A 1 3 ,Ã 1 3 ,Â 1 3,f , andÂ 1 3,F were determined in [24]. The unintegrated tree-level massive quark-antiquark antenna A 0 3 was already introduced in section 2.2. The Born times one-loop interference term δM 1 2 is defined in (2.9).
Compensation term for oversubtracted poles. In certain regions of phase space, the subtraction terms (3.20) and (3.22) exhibit IR singularities that do not coincide with respective singularities in the unsubtracted real-virtual cross section (3.18). In order to remove these spurious singularities one has to introduce an additional subtraction term that is given by [24] dσ T,c,QQg (3.23) The integrands are given in (3.9) and (3.17), respectively. The integration over the phase space of one unresolved parton yields

JHEP12(2016)098
Summary. Combining eqs. (3.18), (3.20), (3.22), and (3.24) yields an expression that is free of (explicit and implicit) singularities in the entire three-parton phase space in D = 4 dimensions: We recall that the terms dσ T,a,QQg NNLO and dσ T,c,QQg NNLO are counterbalanced by the double-real subtraction terms dσ S,a,QQij NNLO and dσ S,b,1,QQij NNLO (ij = gg, qq), respectively, that were defined in section 3.1. Hence, only the integrated form of dσ T,b,QQg NNLO has to be added back to the double virtual contribution that will be discussed in the next subsection.

Double virtual corrections
Finally, we discuss how to compute the order α 2 s contribution of the QQ final state to the differential massive quark-pair production cross section within the antenna framework, that is, the sum of the three terms in the second line of (3.1).
Unsubtracted real-virtual cross section. The renormalized one-loop and 2-loop QQ matrix elements defined in (2.10) yield the following O(α 2 s ) correction to the differential cross section: Summation over all colors and spins is understood. The spin averaging factor 1/4 for unpolarized e − e + collisions is included in (3.26). The factorC( ) is defined in (2.8).
The real and imaginary parts of the one-loop vertex functions V QQ (V = γ, Z), up to and including terms of order , and the corresponding two-loop vertex functions to order 0 were computed in [51][52][53]; cf. also [54]. With these vertex functions, (3.26) can be computed in straightforward fashion. The last term on the right-hand side of (3.26) can be decomposed into different color structures, that is, leading and subleading color contributions, terms that involve a massless and massive quark loop in the gluon vacuum polarization tensor, and triangle contributions summed over all quark flavors. These triangle contributions, which are finite [53], are non-universal QCD corrections to the leading-order cross section.
Subtraction term. Recalling the subtraction terms that were introduced above, those that remain to be counterbalanced are dσ T,b,QQg NNLO (cf. eq. (3.22)) and dσ S,b,2,QQij (ij = qq, gg) (cf. eq. (3.8) and (3.16)). They have to be integrated over the unresolved one-parton, respectively two-parton phase space in order to serve as counterterm for subtracting the

JHEP12(2016)098
IR poles in of the double-virtual correction (3.26). We get The variable y is defined in (2.14). The integrated antenna functions B 0 4 , A 0 4 , andÃ 0 4 were computed in [38,39] and A 1 3 ,Ã 1 3 ,Â 1 3,f , andÂ 1 3,F were determined in [24]. The subtraction term (3.27) has to be added to (3.26). In the sum all IR poles cancel in analytic fashion. After summing these terms and after analytic cancellation of the IR poles, one can take the limit → 0 and perform the remaining integration over the two-parton phase space in four dimensions.

Results for top-quark pair production
In this section we present our numerical results for the total tt cross section and for several distributions, including the top-quark forward-backward asymmetry above the tt threshold at order α 2 s . We use the input values m W = 80.385 GeV, m Z = 91.1876 GeV, and Γ Z = 2.4952 GeV [55]. We use m t = 173.34 GeV for the top-quark mass in the on-shell scheme. The other quarks are taken to be massless. The sine of the weak mixing angle, s W , is fixed by For computing the electroweak couplings we use the socalled G µ scheme (cf., for instance, [14]) where the electromagnetic coupling is given by α = √ 2G µ m 2 W s 2 W /π = 7.5624 × 10 −3 with G µ = 1.166379 × 10 −5 GeV −2 . The running of the MS QCD coupling α s (µ) is determined in f = 6 flavor QCD from the input value α (f =5) s (µ = m Z ) = 0.118. In this section µ refers to the renormalization scale. Because we work to lowest order in the electroweak couplings each of the various contributions dσ (i,j) to the differential QQ cross section at order α 2 s discussed in section 2

JHEP12(2016)098
and 3 is given by the sum of an s-channel γ and Z-boson contribution and a γZ interference term. The dσ (i,j) have the structure Here L µν a denote the lepton tensors (with the boson propagators included) and H a is necessary because some second-order matrix elements involve also terms that contain the electroweak couplings of q = Q, see below.) In this work we consider unpolarized e − e + collisions. We separate each contribution (i, j) on the right-hand side of (4.1) into a parity-even and -odd term.
The two-loop tt matrix elements and the integrated antenna subtraction terms discussed in section 3 contain harmonic polylogarithms (HPL) [56]. We evaluate them with the codes of refs. [57,58]. The integrated antenna functions A 1 3 ,Ã 1 3 that appear in (3.27) are expressed in terms of HPL and cyclotomic harmonic polylogarithms [59][60][61]. We evaluate them numerically by using the integral representation of these functions.
For center-of-mass (c.m.) energies √ s > 4m t , four-top production, i.e., tttt production occurs. The order α 2 s cross section of this process is infrared-finite. It makes only a small contribution to the inclusive tt cross section. Moreover, the tttt final state has a distinct signature and could be experimentally distinguished from tt final states. Below we consider c.m. energies √ s 4m t .

Cross section and distributions
The order α s and α 2 s corrections ∆ 1 and ∆ 2 are displayed, for three scales µ, in figure 2 as functions of the c.m. energy √ s. The changes of ∆ 1 and ∆ 2 due to scale variations are small -notice the logarithmic scale on the y axis of figure 2.
We have included in the computation of σ NNLO also the non-universal contributions of order O(α 2 s ) (cf. section 3) that contain the electroweak couplings of quarks q = t. These contributions are, however, very small. For instance, at √ s = 500 GeV they amount to −0.16% of the total second order correction ∆ 2 defined in (4.2), and this fraction decreases in magnitude for smaller c.m. energies. Our results displayed in figure 2 agree with the calculation of the tt cross section in [26], shown for µ = √ s in figure 1 of this reference. Moreover, considering only the electroweak vector-current contributions to σ NNLO , we agree with the results of [24,25]. In addition, we have compared also with the analytically known threshold expansions [62][63][64][65] and asymptotic expansions [19][20][21][22] of σ NNLO in the regimes α s β 1 (where β is the top-quark velocity) and m 2 t /s 1, respectively, and find agreement. Close to the tt threshold the fixed order perturbative expansion of the cross section and distributions breaks down due to Coulomb singularities. (This kinematic regime has been analyzed in detail with effective field methods.) One can see in figure 1 the onset of the 1/β singularity in the NNLO cross section for √ s → 2m t . We list in table 1 the QCD corrections ∆ 1 and ∆ 2 for selected c.m. energies √ s for µ = √ s. With the input values given above the cross section σ NNLO reaches its maximum at √ s = 381.3 GeV. We obtain σ NNLO (381.3GeV) = 0.843 pb for µ = √ s. The numbers in table 1 and figure 2 suggest that fixed order perturbation theory can be applied for √ s > 360 GeV. Next we turn to differential distributions. We consider the distribution of the cosine of the top-quark scattering angle θ t = ∠(t, e − ) in the c.m. frame, the transverse momentum p t T of the top quark and of the tt system, p tt T = |k T,t +k T,t | with respect to the beam direction, and of the tt invariant mass distribution M tt . In the following we use the schematic notation LO, NLO, and NNLO for dσ LO /dO, dσ NLO /dO = (dσ LO + dσ 1 )/dO, and dσ NNLO /dO = JHEP12(2016)098  (dσ LO + dσ 1 + dσ 2 )/dO, where O denotes one of these observables. We confine ourselves to c.m. energies 400 and 500 GeV where the tt cross section is rather large.

JHEP12(2016)098
The plots in figure 3 display the distribution of cos θ t at √ s = 400 and 500 GeV at LO, NLO, and NNLO QCD. In the upper plots, which refer to µ = √ s, the NLO and NNLO distribution is given by the sum of the grey and red band and by the sum of the grey, red, and blue band, respectively. As expected the first-and second-order QCD corrections decrease if one moves further away from threshold. As the panels in the middle of the plots show the inclusion of the order α 2 s correction significantly reduces the dependence of the distribution on variations of the scale. Both the order α s and order α 2 s follow the same pattern as the leading-order distribution: they are larger in the top-quark forward direction and thus increase the top-quark forward-backward asymmetry, cf. section 4.2. The ratios dσ 1 /dσ LO and dσ 2 /dσ LO for √ s = 400 GeV and µ = √ s shown in the lower panel of the left plot in figure 3 agree with the corresponding result given in [26].   figure 4 agree with the corresponding plot displayed in [26], except for the last bin at (p t T ) max . This is apparently due to the fact that a slightly lower value of the top-quark mass was used in [26]. This shifts (p t T ) max at 400 GeV to a value slightly above 100 GeV.

JHEP12(2016)098
The left plots of figure 5 show, for √ s = 500 GeV, the distribution of the transverse momentum of the tt system, p tt T , for events with p tt T ≥ 10 GeV. The p tt T cut removes the LO QCD contribution and events with very soft massless parton radiation at order α s and α 2 s . For √ s = 500 GeV the maximum p tt T is 129.81 GeV, but events with p tt T near this value are very rare. The NLO and NNLO QCD corrections increase significantly towards small p tt T . This is due to logarithmic enhancement in the variable p tt T that arises in the sum of the order α 2 s three-parton and four-parton contributions. In the bin 10 GeV ≤ p tt T ≤ 20 GeV the order α 2 s correction is almost 50% of the NLO correction. The fixed-order calculation of the distribution becomes unreliable for small p tt T ; the logarithms should be resummed. But this is beyond the scope of this paper. An analogous statement applies to the right plots of figure 5 that show the tt invariant mass distribution for events with M tt ≤ 490 GeV. This cut removes the LO QCD contribution and events with very soft parton radiation.

The forward-backward asymmetry
The top-quark forward-backward asymmetry A FB is defined as the number of t quarks observed in the forward hemisphere minus the number of t quarks in the backward hemisphere, divided by the total number of observed t quarks. Forward and backward hemispheres are defined with respect to a certain IR-safe axis. For top-quark production in e − e + collisions, the top-quark direction of flight is a good choice, because this direction can be reconstructed for instance with lepton plus jets events from tt decay. This axis is infrared-and collinear-safe for massive quarks. Thus A FB is computable in perturbation theory. As long as we consider A FB below the four-top threshold, i.e., as long as we do not include the tttt final state in the computation of the forward-backward asymmetry, A FB can be expressed in terms of the antisymmetric and symmetric tt cross section where the antisymmetric cross section σ A is defined by

JHEP12(2016)098
and σ S is equal to the cross section calculated in section 4.1. As in the previous subsection θ t denotes the angle between the incoming electron and the top-quark direction of flight in the e − e + c.m. frame. The NLO QCD and electroweak corrections to A FB for massive quarks were determined in [66][67][68] and [14,69], respectively. The order α 2 s corrections to A FB were calculated in the limit of massless quarks in [70][71][72][73] and for top quarks with full mass dependence in [26,74].
In order to compute (4.4) to second order in α s in the spirit of perturbation theory, we Taylor-expand (4.4) to second order in α s and obtain where A LO FB is the forward-backward asymmetry at Born level, and A 1 and A 2 are the QCD corrections of O(α s ) and O(α 2 s ), respectively. These terms are given by In analogy to the notation in (4.1) the first number i in the superscript (i, j) labels the number of partons in the final state and the second one the order of α s . Table 2 contains our results for the top-quark forward-backward asymmetry using the expansion (4.6)-(4.9) for several c.m. energies and for the input values as given at the beginning of section 4. Notice that A NLO FB = A LO FB (1 + A 1 ). The central values refer to the scale choice µ = √ s and the given uncertainties are obtained by varying µ between √ s/2 and 2 √ s. We have included in A 2 also the non-universal contributions that contain the electroweak couplings of quarks q = t. (We remark that the square of the diagrams where γ * /Z * couple to q = t and the tt pair is produced by gluon splitting does not contribute to the antisymmetric cross section.) These contributions are, however, small. The ratio r = A non 2 /A 2 of the non-universal and the total order α 2 s correction increases with increasing c.m. energy. We have r = −0.16%, −1%, and −2.4% for √ s = 400 GeV, 500 GeV, and 700 GeV, respectively.
As one can see from table 2, close to the tt threshold, at √ s 360 GeV, fixed order perturbation theory is no longer reliable because the second order correction A 2 to the forwardbackward asymmetry is larger than the first order correction A 1 . For √ s > 380 GeV, the ratio |A 2 /A 1 | becomes smaller than one. Notice that the order α 2 s correction is significant as compared to the first order one: |A 2 /A 1 | 0.5 for the c.m. energies listed in A NNLO FB increases if the top-quark mass is decreased and vice versa. One expects that the top-quark mass can be measured with a much smaller uncertainty than ±0.5 GeV from a tt threshold scan at a future e − e + collider [4,5].
The two-parton, i.e., the tt contribution to A FB is separately IR-finite, both at order α s and at order α 2 s [65]. In the range of c.m. energies given in table 2, the tt final state makes the largest contribution both to A 1 and A 2 . For √ s 500 GeV it is significantly larger than the respective contribution from the three-and four-parton final states. Here, we have computed the tt contribution to A 1 and A 2 with the antenna-subtracted two-parton matrix elements of sections 2 and 3, while it was computed in [65] with the unsubtracted tt matrix elements. We agree with the results of [65]. This serves as a check of our calculation. The sum of the three-and four-parton contribution A 2 to A 2 is also IR-finite. It was computed in [74] with an NLO subtraction scheme, namely dipole subtraction with massive quarks [75]. We agree with the results of [74].
Measurements of forward-backward asymmetries or simulations with Monte Carlo event generators correspond to computations where the ratio in (4.4) is not Taylorexpanded. Using our results for σ A and σ S at O(α s ) and O(α 2 s ), respectively, and the input values as given at the beginning of section 4, we give in table 3 the values for the unexpanded version of the forward-backward asymmetry at NLO and NNLO QCD. (For ease of notation, we use the same symbols as in table 2.) The central values and the uncertainties refer again to the scales µ = √ s and µ = √ s/2 and 2 √ s, respectively. We vary µ simultaneously in the numerator and denominator of (4.4).
The top-quark forward-backward asymmetry at NNLO QCD was computed before in ref. [26] in the unexpanded version with values of m t and α s that differ slightly from the ones that we use here. Our results of table 3 agree with those given in table 1 of that reference.
One may take the spread between the values of the expanded and unexpanded versions of A NLO F B and A NNLO F B given in tables 2, 3 as an estimate of the uncalculated higher order corrections. Alternatively, one may add the scale uncertainties and the uncertainties due to δm t = ±0.5 GeV of the expanded version (cf.  cantly smaller than the projected experimental precision of top-quark A FB measurements at future electron-positron colliders [77,81]. This observable has a high sensitivity to precisely determine the neutral current couplings of the top quark and probe for anomalous couplings [77][78][79][80][81].

Summary
We have formulated, within the antenna subtraction framework, the set-up for calculating the production of a massive quark-antiquark pair in electron-positron collisions at NNLO QCD. Our approach is fully differential in the phase-space variables and can be used to compute any infrared-safe observable. We have applied this formalism to tt production in the continuum and we have calculated, besides the tt cross section also several distributions in order to signify the usefulness of this approach, namely the cos θ t and transverse momentum distribution of the top quark, the transverse momentum of the tt system and the tt invariant mass distribution. The NNLO QCD corrections are sizable for c.m. energies not too far away from the tt threshold. We have also computed the top-quark forwardbackward asymmetry, which is an important observable for determining the neutral-current couplings of the top quark at future lepton colliders, at NNLO QCD. Our result agrees with previous calculations [26,74] of this asymmetry at order α 2 s . Our set-up may be used to investigate a number of other reactions at NNLO QCD where a massive quark-pair is produced by an uncolored initial state. Of interest for future lepton colliders would be the production of tt pairs with spin correlations included. Other applications include the production of charm and bottom quarks, in particular at the Zboson resonance.

A Phase-space mappings that involve massive particles
We describe here the phase-space mappings that are used in the construction of the antenna subtraction terms of section 2 and 3. The momentum mappings required in our case are related either to a single or double unresolved parton configuration in the final state. These mappings must obey four-momentum conservation, must keep the mapped momenta on their respective mass shell, and the mapped momenta must converge to the correct momentum configurations in the the soft and collinear limits. We follow the mapping JHEP12(2016)098 procedures of [48], which apply to the case where all partons are massless, respectively of [33,46] where the massless case was extended to configurations involving massive partons. The analytic formulas of Abelof and Gehrmann-De Ridder that keep the mapped momenta on-shell in the massive case have not been published so far [46]. Therefore we describe below an alternative mapping method for computing the observables used in this paper.

A.1 Three parton final states
We consider the final state Q(k 1 )Q(k 2 )g(k 3 ). The NLO subtraction term dσ S QQg of eq. (2.12) and the NNLO subtraction terms dσ T,b,QQg NNLO and dσ T,c,QQg NNLO of eq. (3.22) and (3.24) depend on mapped momenta obtained from a 3 → 2 mapping k 1 , k 3 , k 2 → k 13 , k 32 . Let's consider the mapping k 1 , k 3 , k 2 → p I ≡ p 13 , p J ≡ p 32 defined in [48] and in appendix B1.1 of [33]: where the parameters x, r, z are given in [33,48]. The mapping (A.1) satisfies fourmomentum conservation, p I + p J = k 1 + k 3 + k 2 , and the mapped momenta behave correctly when the gluon becomes soft: p I → k 1 , p J → k 2 if k 3 → 0. If all three partons were massless, the mapped momenta remain massless, p 2 I = p 2 J = 0 [48]. However, for a massive quark Q modified formulas must be used for the parameters x, r, z in order to get p 2 I , p 2 J = m 2 Q [46]. We recall that on-shellness of the mapped momenta is crucial for deriving the correct integrated antenna subtraction terms from the unintegrated ones.
Here we describe, as an alternative to the analytic formulas of Abelof and Gehrmann-De Ridder [46], a numerical method to obtain on-shell mapped momenta k 13 , k 32 . We use the mapping (A.1) with the parameters x, r, z given in [33] for an intermediate step.
Four-momentum conservation in the e − e + c.m. frame reads: The second equation is the crucial one. It allows to rescale the 3-momenta by a factor ξ such that the 4-momenta p µ I , p µ J are transformed into on-shell 4-momenta k I , k J with mass m I = m J = m Q without destroying 4-momentum conservation.

JHEP12(2016)098
For the second set of momenta in (A.5) mapped on-shell momenta are constructed in completely analogous fashion with the 'spectator' k 2 replaced by k 1 . The above procedure applies also to the various 3 → 2 mappings that are required for the a-type subtraction term dσ S,a,QQgg NNLO of (3.15). Abelof and Gehrmann-De Ridder have derived analytic formulas for the mapped onshell 4-momenta k µ I , k µ J [46]. 4 → 2 mappings for b-type antenna subtraction terms. The b-type antenna subtraction term dσ S,b,2,QQqq NNLO of (3.8) is evaluated with mapped momenta that are obtained by a 4 → 2 mapping k 1 , k 3 , k 4 , k 2 → k I ≡ k 134 , k J ≡ k 342 . (A.11) As an intermediate step we use the 4 → 2 mapping k 1 , k 3 , k 4 , k 2 → p I , p J , where 12) and the parameters x, r 1 , r 2 , z are given 1 in appendix B.2.1 of [33]. For general configurations k j the mapped momenta are not on the m Q mass shell, p 2 I , p 2 J = m 2 Q . Four-momentum conservation in the e − e + c.m. frame reads √ s = p 0I + p 0J , 0 = p I + p J . Let's consider (A.16). The 3 → 2 mapping I is done as described below eq. (A.5): boost to the rest frame of Q µ 2 , rescale, and then boost back to the e − e + c.m. frame. The rescaling involved in the subsequent mapping II is done directly in the e − e + c.m. frame. This yields the two mapped on-shell momenta with mass m Q on the right-hand side of (A.16). The iterated 3 → 2 mappings (A.17) and those involved in constructing the antenna subtraction term dσ S,b,1,QQgg NNLO of (3.17) are performed in analogous fashion.