The next-to-next-to-leading order soft function for top quark pair production

We present the first calculation of the next-to-next-to-leading order threshold soft function for top quark pair production at hadron colliders, with full velocity dependence of the massive top quarks. Our results are fully analytic, and can be entirely written in terms of generalized polylogarithms. The scale-dependence of our result coincides with the well-known two-loop anomalous dimension matrix including the three-parton correlations, which at the two-loop order only appear when more than one massive partons are involved in the scattering process. In the boosted limit, our result exhibits the expected factorization property of mass logarithms, which leads to a consistent extraction of the soft fragmentation function. The next-to-next-to-leading order soft function obtained in this paper is an important ingredient for threshold resummation at the next-to-next-to-next-to-leading logarithmic accuracy.


Introduction
Top quark pair production is one of the most important processes at the Large Hadron Collider (LHC). The total cross section at √ s = 13 TeV is about 800 pb. Both the ATLAS and the CMS experiments have collected nearly 100 fb −1 of integrated luminosity at 13 TeV. This corresponds to 160 million tt events in total. With such a large event sample, top quark physics has become one of the precision frontier of particle physics. Many important measurements related to the top quark, e.g., inclusive and differential cross sections [1,2], top-quark mass and width [3,4], top-quark polarization [5] and et al., can now be done at unprecedented precisions. The large production cross section also allows precision measurements of boosted top quark pairs, which is important for high-mass tt resonances search [6], and for precision study of boosted top-quark jet [7]. Currently, the best fixed-order calculations for top-quark pair production is at the nextto-next-to-leading order (NNLO) in QCD [8][9][10][11][12][13][14] and the next-to-leading order (NLO) in the electroweak coupling [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Recently, these two corrections have been combined to give a more comprehensive description of tt production in [30]. Despite the high precisions of these perturbative calculations, the complicated kinematics of tt production makes it necessary to consider even higher order corrections. This is particularly important since the high energy of the LHC has opened up the possibility to produce "boosted" top quark pairs, which means that the energies of the top quarks are much larger than their rest mass m t . In [14], it has been found that the NNLO QCD differential cross sections in the boosted regime are rather sensitive to the choice of factorization and renormalization scales. This scale dependence can be dramatically reduced by resumming certain towers of large logarithms to all orders in the strong coupling α s [31,32]. These include not only the threshold logarithms which arise when the partonic center-of-mass energy approaches the tt invariant mass M, but also the mass logarithms of the form ln n (m 2 t /M 2 ) which develop in the boosted region M ≫ m t . The resummation of the threshold logarithms in tt production requires a couple of ingredients, such as the hard function and the soft function, as well as various anomalous dimensions. The NLO hard and soft functions have been computed in [33], and the anomalous dimension matrices have been derived in [34,35]. These enabled the resummation to be done at the nextto-next-to-leading logarithmic (NNLL) accuracy [33]. Given the NNLO accuracy achieved from fixed-order calculations, it is desirable to extend the threshold resummation to N 3 LL. Such a calculation would improve the theoretical predictions over the whole phase space, all the way from low invariant mass to the boost regime. In order to achieve that, the NNLO hard and soft functions are necessary. The NNLO hard function can in principle be extracted from the virtual amplitude calculated in [36]. Therefore, the NNLO soft function becomes a major bottleneck in pushing up the resummation accuracy of tt production, and is the subject of this article.
The soft functions describe the cross sections in the soft limit. The behavior of scattering amplitudes and cross sections in the soft limit is of high interest not only phenomenologically, but also theoretically. For example, the soft theorems in gauge theories and in gravitational theories [37][38][39] are of fundamental importance in the understanding of their structures. In perturbative calculations in gauge theories, both the exchange of virtual soft particles and the emission of real soft ones can lead to infrared (IR) divergences. These must cancel against each other in order to arrive at meaningful predictions for physical observables. While such cancellations have been proven generically [40,41], the practical treatments of the IR divergences are highly non-trivial. Both the virtual and real contributions need to be calculated analytically in order to verify the precise cancellation. For the virtual amplitudes, when all external hard partons are massless, the soft singularities enjoy a dipole form up to two loops [42], thanks to the emerged scaling symmetry as the energy of hard partons become large [43,44]. Non-trivial corrections to the dipole form of soft singularities for massless scattering first appear at three loops, and have been computed recently in [45]. The situation for massive amplitudes is much more complicated. Non-trivial correlations among three or four partons appear already at two loops [34,35,[46][47][48]. These virtual singularities must have the same structure as the real ones, and the soft functions provide a perfect place to investigate the latter. It is therefore interesting to calculate the massive soft functions through to NNLO, in order to study these multi-parton correlations from a different perspective.
The soft functions are defined as the vacuum expectation values of certain operators consisting of light-like and time-like soft Wilson lines. In simpler situations, they have been extensively studied in the literature. For processes involving two massless partons, such as the Drell-Yan process and Higgs production through gluon fusion, the soft functions have been calculated up to the N 3 LO [49]. For processes with 4 massless partons such as di-jet production and boosted heavy quark pair production, the NNLO soft function was obtained in [50]. When massive partons are involved, the calculation gets much more complicated. The soft function for the e + e − → tt process has been calculated at the NNLO in [51]. Much less is known in the case of hadronic production of top quark pairs, for which only the NLO soft function is available [33]. Our result in this work therefore serves as the first example of an NNLO soft function for massive scattering with 4 external partons.
This paper is organized as follows. In Section 2 we lay out the generic definition and renormalization of threshold soft functions. In Section 3 we provide the result of the NLO soft function to higher powers of the dimensional regulator ǫ, which is a necessary ingredient in the renormalization of the soft function at NNLO. Section 4 describes the main efforts of this work, namely the calculation of the NNLO bare soft function. We then perform its renormalization in Section 5, and discuss some cross-checks and the numerical impact of our new result. Finally, we conclude and discuss some future applications and extensions of our calculation in Section 6.

Formalism
We consider a generic scattering process involving energetic massless quarks, gluons and massive partons (such as top quarks or some new colored particles often present in models beyond the SM). The interactions of soft gluons with these energetic partons can be described by Wilson lines defined as where P denotes path ordering, v i is a 4-vector pointing to the direction of the momentum of the i-th parton, which satisfies v 2 i = 0 for massless partons and v 2 i > 0 for massive partons.
Note that here we have taken all vectors v i to be incoming. The boldface T a i is the color generator associated with the i-th parton in the color-space formalism [52,53]. It is evident that the Wilson lines are invariant under the rescaling v i → λv i for any λ > 0, since this change can be compensated by a change of the integration variable s → s/λ. We could employ this freedom to normalize the direction vectors of massive partons to v 2 i = 1. This has the physical meaning that v i is the 4-velocity of the i-th parton: v i = ±p i /m i where m i is the mass of the i-th parton. However, we'd like to keep this possibility open for the sake of generality.
Putting the Wilson lines together, the behavior of the n-parton scattering amplitude in the soft limit can be obtained via studying the vacuum matrix elements of the Wilson loop operator constructed out of the Wilson lines where {v} denotes the collection of the directional vectors v i , x is a time-like vector, and T denotes time-ordering. It is well-known that the vacuum matrix elements of the Wilson loop operator, when calculated in perturbation theory, contain ultraviolet (UV) divergences which need to be renormalized [54,55]. The renormalization properties of the Wilson loops can be used to study the infrared singularities of scattering amplitudes, as was illustrated in [34,35,43,47,56]. The Wilson loop operator is also an essential ingredient in the factorization of scattering cross sections in the soft limit. We consider scattering processes at hadron-hadron colliders with no final state massless partons at the leading order. These include, for example, top quark pair production (possibly associated with other colorless particles such as the Higgs boson and electroweak gauge bosons), production of 4 top quarks, squark and gluino productions in supersymmetric models, as well as productions of top partners in many new physics models. At higher orders in the strong coupling constant, there will be additional emissions of gluons and quarks in the final state. We are interested in the case where these additional emissions are all soft, i.e., with energies much smaller than the typical momentum transfer of the hard-scattering process. Note that the precise meaning of "soft" depends on the reference frame, which leads to different forms of the factorization formula, such as the "pair-invariant-mass" (PIM) kinematics and the "single-particle-inclusive" (1PI) kinematics in top quark pair production discussed, e.g., in [33,57,58]. While the formalism can be applied to any reference frame, in the following, we will work in the center-of-mass frame of the two incoming partons, which are not only good for demonstration purposes, and are also adopted in many existing calculations. For example, this corresponds to the PIM kinematics in [33,59,60] for tt production, in [61,62] for stop pair production, and in [63,64] for ttH production. Schematically, we are considering partonic processes of the form where h 1 and h 2 are two incoming massless partons, h I (I = 3, . . . , n) are outgoing massive partons, X denotes colorless particles such as the Higgs boson and electroweak gauge bosons, and X s represents additional soft radiations which we want to describe. Here, the momenta p I (I = 3, . . . , n) are chosen to be outgoing, but we still keep v I to be incoming for convenience.
In the center-of-mass frame of the two incoming partons, the emissions of additional soft partons are described by the so-called soft function S, which is simply the momentum-space version of the vacuum matrix element of the Wilson loop operator: where the reference vector v 0 = (2, 0, 0, 0), and ω represents (2 times) the energy of the additional soft partons. We have included a normalization factor such that the definition of the soft function coincides with that in [33]. Here d 1 and d 2 are the dimensions of the SU(N) color representations to which the partons h 1 and h 2 belong. For later convenience, it is useful to perform a Mellin or Laplace transform into the moment spacẽ where Λ is a soft momentum scale in the moment space. As discussed above, the bare soft function contains UV divergences which can be regularized in dimensional regularization with d = 4 − 2ǫ. These UV divergences are removed by the operator renormalizatioñ where µ is the renormalization scale and L ≡ ln(Λ 2 /µ 2 ). As indicated in the above formula, both the renormalized soft functions(L, {v}, µ) and the renormalization factor Z s (L, {v}, µ) are µ-dependent and satisfy renormalization group equations (RGEs) The generic form of the soft anomalous dimension matrix Γ s up to the two-loop order can be written as [47] where the lower-case indices (i, j, k) run over massless partons 1 and 2, while the capital indices (I, J, K) run over massive partons. We have introduced the abbreviations L s ≡ L − iπ and w Ij ≡ v I · v j + i0. The notation (I, J, . . .) denotes unordered tuples of distinct parton indices. The functions γ cusp (α s ) and γ cusp (β IJ , α s ) are the famous cusp anomalous dimensions for lightlike Wilson lines and time-like Wilson lines, respectively [54,65,66] γ cusp (α s ) = α s π + α s 4π with the cusp angle given by The single parton soft anomalous dimensions γ i s (α s ) and γ I (α s ) are where C i(I) = C F for the fundamental representation, and C i(I) = C A for the adjoint representation of the gauge group. The three-parton correlation functions F 1 and f 2 were calculated in [34,35]. The function F 1 describes correlations among three massive partons, and can be written as where the indices (L,M,N) run over (I,J,K) with ǫ LM N = 1 if (L,M,N) is an even permutation of (I,J,K), and The function f 2 describes correlations among two massive partons and one massless parton, and is given by Given the anomalous dimension matrix Γ s , one can solve the RGE (7) to obtain the renormalization factor Z s . To this end, it is useful to decompose Γ s in (9) into the form and then 3 The soft function for tt production and the NLO result to arbitrary orders in ǫ While the formalism introduced in the last section is very generic and applies to a lot of processes, the actual calculation of the soft function could get very complicated when the number of independent scalar products v i · v j becomes large. In this paper, we begin with the special case of tt production, where the partonic processes can be described as where p 2 3 = p 2 4 = m 2 t with m t the mass of the top quark. In the soft limit p s → 0, there are 3 independent Lorentz invariant kinematic variables, which can be chosen as m t and It is convenient to introduce dimensionless quantities β and y = cos θ, defined as where β and θ have the physical meanings of the 3-velocity and the scattering angle of the top quark in the partonic center-of-mass frame. The soft function depends on β and y through the following combinations of the directional vectors v i : and we also have To calculate the bare soft function, it is convenient to start from the momentum-space version S bare (ω, β, y) introduced in Eq. (4). Here we have expressed the dependence on the directional vectors v i through the quantities β and y. The perturbative expansion of the momentum-space soft function can be written as where 1 denotes the identity operator in color space, and the NLO soft function S bare (ω, β, y) was already calculated in [33]. In fact, since the NLO soft function involves 2-parton correlations at most, the same calculation can be applied to scattering processes with more than two colored particles in the final state. This has been done, e.g., for the case of ttH production in [63,64], and can be extended to more complicated processes such the simultaneous production of two tt pairs. In order to calculate the NNLO soft function, however, we will need the NLO one to higher orders in the dimensional regulator ǫ, which will produce a finite contribution to the renormalized NNLO soft function. We will describe such a calculation in the following, and the calculation of the NNLO bare soft function will be discussed in Section 4.
In [33], the NLO soft function was calculated by brute-force evaluation of the relevant phase-space integrals. While it is possible to continue using such a method to obtain the higher order terms in ǫ, it is useful to employ a more systematic approach which can be extended to the NNLO calculation.
The definition (4) of the soft function involves a summation over soft final states X s . It is easy to see that when |X s is the vacuum state, i.e., when there is no extra soft emission, the matrix elements involve scaleless integrals in dimensional regularization, which are defined to be zero. Therefore, at the NLO, the only contribution is given by where a summation over the helicity and the color of the gluon is understood. We use QGRAF [67] to generate the squared amplitudes in the above formula, and use FORM [68] to manipulate the resulting expression. The phase-space integrals appearing in the result have the general form . From symmetry considerations, it is obvious that we only need to calculate I 12 , I 13 , I 33 and I 34 . At this point, we note that while the soft function itself does not depend on the normalizations of the directional vectors v i , it is convenient to fix them in practical calculations. Therefore we will choose the normalizations v and v 2 3 = v 2 4 = 1 − β 2 in the following. Note that the normalizations of v 3 and v 4 are unconventional. In the center-of-mass frame of the incoming partons, these vectors are then parameterized by The phase-space integrals appearing in the result are not independent, and we employ the integration-by-parts (IBP) [69,70] method to find relations among them. To this end it is necessary to use the relation which is known as the reverse unitarity method [71], to express the phase-space integrals in terms of loop integrals. The δ-function in Eq. (27) can be similarly written as These integrals are then feed into the program packages Reduze2 [72] and FIRE5 [73], which use the IBP relations to reduce the relevant loop integrals to a number of master integrals. After the IBP reduction, one can recover the phase-space integrals by reversing the relations (29) and (30). In the NLO case, the master integrals can be chosen as F 0,0 , F 0,1 and F 1,1 , where F a 1 ,a 2 is defined as In order to calculate the master integrals to arbitrary orders in the dimensional regulator ǫ, we employ the method of differential equations [74,75]. Taking the partial derivative of F a 1 ,a 2 with respect to β will lead to integrals with a 2 index shifted. However, all these integrals can be expressed as linear combinations of the master integrals. We collect the master integrals into a vector with 3 components The partial derivative of f with respect to β then has the form whereÂ is a 3 × 3 matrix. The matrixÂ is a rather complicated function of ǫ, β, y and ω, which makes the differential equation not so straightforward to solve. It is possible to simplify the above equation by a linear transformation g(ǫ, β, y) =T (ǫ, β, y, ω) f(ǫ, β, y, ω), where the matrixT is given bŷ The new vector g (so-called "canonical basis") satisfies a simpler differential equation [76] ∂ β g(ǫ, β, y) = ǫB(β, y) g(ǫ, β, y) .
Note that the matrixB does not depend on ǫ and ω anymore, and is given bŷ The vector g has no singularity in ǫ and can be Taylor-expanded in the form The differential equation (35) can then be solved order by order in ǫ: where the boundary conditions g (n) (β 0 , y) at some boundary point β 0 need to be obtained through other methods, which will be discussed later. Such iterated integrals give rise to so-called generalized polylogarithms (GPLs) [77], which are defined by with G(; β) ≡ 1. The special case where all the a i 's are zero is defined as The GPLs have many good mathematical properties (for a review, see e.g. [78]), and can be straightforwardly evaluated by program packages such as GiNaC [79]. They therefore form a wonderful basis for expressing our results.
In order to solve the differential equations (35), we also need the explicit expression of g(ǫ, β, y) at the boundary point β 0 , serving as the boundary condition. In our case, it is convenient to choose the point β = 0 as the boundary. The calculation of the boundary condition can be simplified by observing that the matrixB contains a singular term proportional to 1/β. It is clear that ∂ β F a 1 ,a 2 can produce a 1/β coefficient only if F a 1 ,a 2 itself develops a power-like or logarithmic divergence when β → 0. One can easily check that all the master integrals in f are regular in the limit β → 0. The same applies to the components of g since the transformation matrixT is also regular. It follows that which leads to the conditions with g i being the i-th component of g. We now only need to directly evaluate the component g 1 (which actually doesn't depend on β) at the boundary β = 0, which is very simple: We now have everything we need to express the NLO soft function as an abstract matrix in color space, in terms of the inner products of color generators T i · T j ≡ T a i T a j . This abstract form is generic and especially useful if we want to apply our result to more complicated processes. However, for practical computations of tt cross sections, it is convenient to choose a color basis and express the soft function as a 2 × 2 matrix in the quark-anti-quark annihilation channel, and a 3 × 3 matrix in the gluon fusion channel. Such matrix elements are defined as where {|c I } is an orthogonal color basis. In accordance with [33], we choose the singlet-octet basis with

Double-real contributions
We first present the calculation of the double-real contributions. As in the NLO calculation, we generate relevant Feynman diagrams and amplitudes using QGRAF [67]. The phase-space integrals in the double-real contribution have the generic form where F denotes the integrand consisting of scalar products among the directional vectors v i and the two momenta k 1 and k 2 . We generically call these scalar products "propagators". There exist many different propagators in our squared amplitudes. However, only a subset of them appears in any individual integral. It is therefore useful to classify all the integrals into a couple of "integral families", each defined by a particular set of propagators. For this purpose, we first classify the relevant Feynman diagrams into three categories according to the number of independent Wilson lines involved: 1) those involving one or two Wilson lines; 2) those involving three Wilson lines; and 3) those involving all four Wilson lines. We discuss the calculation of the first two categories in the following. The diagrams involving all four Wilson lines can be trivially expressed as a convolution of two NLO integrals, and we do not bother to discuss them here.
We now turn to the second integral family in the one-or two-Wilson-line diagrams. It is defined by the set of propagators We denote the corresponding integrals as F a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 , defined similar to Eq. (49), but with the propagators D i chosen from the above set (53). The master integrals in this family can be chosen as 0,0,0,0,0,0 , F 1,0,0,0,−1,2 , F 0,0,0,0,1,0 , F 0,0,0,0,1,1 , F 0,0,0,1,0,0 , F 0,0,0,1,1,0 , F 1,0,0,1,1,0 , F 0,0,0,1,1,1 , F which are transformed into a canonical basis g (2) (ǫ, β, y) by the following transformation matrix The differential equation satisfied by g (2) is given by withâ (2) ,b (2) andĉ (2) being 9 × 9 constant matrices. In order to solve the differential equations (52) and (56), we also need the boundary conditions at some value of β. Similar to the NLO case, we again choose the point β = 0 as the boundary, where only 7 of the integrals in g (1) and g (2) are non-vanishing. Some of the boundary conditions are related to each other, similar to (43). The independent ones are given by The hypergeometric function appearing in the above formula can be expanded in ǫ with the help of the program package HypExp [80]. The differential equations of g (1) and g (2) can then be solved order-by-order in ǫ in terms of the GPLs, similar to the NLO case (39).

Three-Wilson-line diagrams
We now turn to the three-Wilson-line diagrams. In the calculation of the two-loop anomalous dimensions and infrared singularities [34,35], it was found that the three-Wilson-line diagrams are the most complicated ones. They give rise to the three-parton functions F 1 and f 2 in Eq. (9). It is therefore highly interesting to see how these complications appear in the calculation of the NNLO soft function. As will be clear below, the genuine three-parton correlations only arise from the virtual-real contributions, and the double-real three-parton contributions can always be expressed as convolutions of two NLO integrals. This can be understood since the soft function is a Hermitian matrix from its definition (4). On the other hand, the genuine three-parton contributions (such as the functions F 1 and f 2 ) multiply the anti-Hermitian color factor if abc T a i T b j T c k . Therefore the soft function can only receive contributions from the imaginary parts of the three-parton integrals, which are only present in the virtual-real diagrams but not in the double-real diagrams.  More practically, one can carry out an analysis similar to that in [50] (done for the massless soft function). In that paper, it was demonstrated that the three-Wilson-line integrals in the double-real contributions can be combined into convolutions of NLO integrals. The same applies to the massive soft function. To see this, consider for example the set of diagrams shown in Figure 3. Each of these diagrams gives rise to rather complicated integrals. However, the sum of the four diagrams leads to the simple integral which is obviously a convolution of two NLO integrals. This fact does not depend on whether v i , v j and v k are light-like or time-like, and therefore applies equally well to the massless and the massive soft functions. Another example is the two diagrams involving the three-gluon vertex, shown in Figure 4. As demonstrated in [50], they add up to zero due to the color structure, irrespective of the nature of the Wilson lines involved.
The 15 × 15 matricesT ′(3) ,â (3) ,b (3) ,ĉ (3) ,d (3) ,ê (3) as well as the 21 × 21 matricesT ′(4) ,â (4) , b (4) ,ĉ (4) ,d (4) ,ê (4) are non-diagonal with lengthy expressions. We choose to give them in an electronic file attached to the arXiv submission, together with the matricesâ (1) ,b (1) ,ĉ (1) ,d (1) , e (1) ,â (2) ,b (2) andĉ (2) appearing in the calculation of the double-real contributions. In order to solve the differential equations (65) and (69), we again need to calculate the boundary conditions at β = 0. In the virtual-real case, an additional complication arises due to the fact that some of the master integrals are logarithmic divergent in the limit β → 0. This happens when the loop integral (involving both v 3 and v 4 ) give rise to Coulomb-like singularities of 1/β which are multiplied by some ǫ-dependent powers of β. When this is the case, we cannot set β = 0 before performing the integration, like what we did in the doublereal case. Instead, we have to explicitly calculate the β-dependence of such master integrals (at least their asymptotic form as β → 0). Fortunately, this only needs to be done for two of the master integrals, which are g (4) 6 (ǫ, β, y) = Here, we explicitly write the imaginary parts of the propagators involving the loop momentum l, since the loop integrals over l depend on them. In the following, we will suppress the imaginary parts and the +i0 prescription is always understood. The results of the integrals over the loop momentum l can be found in [81,82]. The remaining integrations over the momentum of the real gluon k can be carried out in the limit β → 0, since no new divergence arises in this limit. The asymptotic behaviors of g (4) 6 and g (4) 9 near the boundary can then be obtained and are given by This can be readily expanded in ǫ to arrive at the boundary conditions. The other independent boundary conditions are easier to obtain and are given by With the above results, it is straightforward to derive the virtual-real contributions to the NNLO soft function. It should be emphasized again that the three-Wilson-line diagrams give non-vanishing contributions, which are necessary to give the correct pole structure in accordance with the anomalous dimension (9). On the other hand, these contributions are proportional to the anti-Hermitian color factor if abc T a i T b j T c k , and only enter the off-diagonal entries of the soft function in the singlet-octet basis (46). They therefore do not appear at the level of the NNLO cross section, but will enter the resummation formula which encodes higher order effects beyond the NNLO.

The bare soft function in the moment space
Assembling all the ingredients in the last two subsections, we obtain the NNLO momentumspace bare soft function S (2) bare (ω, β, y). It is written in terms of star-distributions in ω. For later convenience, it is useful to transform the soft function to the moment space using a Laplace transform (5). This can be most easily done by observing that the ω-dependence of the NNLO bare soft function comes from an overall factor ω −1−4ǫ , and A similar transformation rule can be derived for the NLO bare soft function. The resulting moment-space soft function can then be written as a function of Λ: Similar to the momentum-space soft function, we define the matrix elements of the momentspace soft function ass where the color basis is chosen as in Eq. (46). This matrix-valued moment-space soft function will be the main object to be studied in the following.

Validations of the results
While our results of the NNLO soft functions are novel, it is possible to partially validate them by checking their consistency with some known results in the literature. We have performed 3 checks: 1) that they satisfy RGEs according to the anomalous dimension matrices calculated in [34,35]; 2) that they reproduce in the threshold limit the results of [83,84]; 3) that they correctly factorize in the boosted limit according to the factorization formula given in [85]. As discussed in Section 2, the renormalized soft function satisfies a renormalization group equation This property is closely related to the renormalization in Eq. (78). Given that our result is correctly renormalized, we expect that it should naturally satisfies Eq. (85). Nevertheless, we have calculated the left side and the right side of the above equation and indeed find consistency.
In the threshold limit s → 4m 2 t or β → 0, the top quark and the anti-top quark are at rest in the partonic center-of-mass frame. In this case, if the tt pair forms a color-singlet state, the soft gluons cannot probe it, and the situation is no different from the Drell-Yan process or the Higgs boson production process. The corresponding element of the soft function matrix then reduce to those calculated in [83]. On the other hand, if the tt pair forms a color-octet state, the corresponding matrix element in the threshold limit has been calculated in [84]. We can therefore check the consistency of our result by taking the limit β → 0 in our expressions. This is easy to do in our formalism since β = 0 is essentially the boundary point of the differential equations. One can also directly take the β → 0 limit starting from the explicit expressions ofs(L, β, cos θ, µ), which is expressed in terms of GPLs of β. Using the property that G(a 1 , . . . , a n ; 0) = 0 unless all the indices are zero, this limit is rather straightforward to obtain. We have done this simple exercise and find that the results are in perfect agreement with those in [83] and [84]. In particular, for the qq channel we find the diagonal entries to bẽ s qq,(2) 11 (L, β → 0, cos θ) The results for the gg channel can be simply obtained from the above expressions by the replacement C F → C A . An interesting subtlety in the above exercise is that the non-diagonal entries of the soft function do not vanish in the limit β → 0. Instead, they develop logarithmic divergent behaviors in that limit. These non-vanishing terms arise from the three-Wilson-line virtualreal diagrams and are therefore purely imaginary. For example, we havẽ s qq,(2) 12 This is actually expected since the f 2 function in the anomalous dimension matrix has the similar property in the threshold limit (see, e.g., Section 3.4 of [35]). However, such a behavior cannot be seen from the calculation of [83,84], and is a novel feature of our results. Besides the above "trivial" checks, a highly non-trivial cross-check of our result is the opposite, boosted limit β → 1 or M 2 ≫ 4m 2 t . In this limit, it was shown in [85] that the soft function should factorize in the form s massive ln Λ 2 µ 2 , β, cos θ, µ =s massless ln Λ 2 µ 2 , cos θ, µ s 2 D ln wheres massive is the massive soft function calculated in this work,s massless is the soft function with top quarks treated as massless (which was calculated at NNLO in [50] 1 ), ands D is a soft-collinear partonic fragmentation function describing a boosted top quark evolving into a top quark plus soft radiations. The soft-collinear fragmentation functions D was not directly calculated in the literature. In [86,87], it was related to the shape function in B-meson decays and was extracted from existing result of the latter. In the Appendix A of [85], however, it was found that this result ofs D is inconsistent with other related calculations. To understand this, we note thats D should satisfy another factorization formulã where L D ≡ ln m 2 t Λ 2 /M 2 µ 2 (note thats D depends on a see-saw scale m t Λ/M). In the above formula,D t/t is the full partonic fragmentation function of the top quark in the moment space, and C D is a hard-collinear matching coefficient. The momentum-space version ofD t/t was calculated in [88]. Using this result and the result ofs D from [86,87], it is possible to extract the form of C D from Eq. (90). However, the function C D could also be extracted from the high-energy behavior of the scattering-amplitude involving heavy quarks [89]. The main conclusion of the Appendix A of [85] is that these two extractions do not coincide! Using our new result of the NNLO massive soft function, it is for the first time that one can directly validate the factorization formula (89) at the NNLO, and extract the soft-collinear fragmentation functions D at this order. It is then possible to resolve the conflict between the two results of C D . In order to do these, we need to take the limit m t → 0 or β → 1, and carefully extract the logarithms of m t . Such logarithms arise from GPLs which are singular in the limit β → 1. In order to extract such singularities, we employ properties of GPLs to convert their argument to cos θ, and put all β dependences into the weights. We have used the program package HyperInt [90] to accomplish this. After this conversion, the singular terms become powers of ln(1 − β) ≈ ln(2m 2 t /M 2 ). Finally, we insert the results and the massless soft function from [50] into Eq. (89), and we find the NNLO coefficient ofs D to bẽ The above formula differs from that in [86,87] by a constant term 4π 2 C A C F , which is essentially the inconsistency discussed in the Appendix A of [85]. Accordingly, the NNLO coefficient of the function C D is given by where L m ≡ ln µ 2 /m 2 t . The above expression of C D satisfies both Eq. (90) and the highenergy limit of the scattering amplitude. The inconsistency on the form of C D found in the Appendix A of [85] is thus resolved by our calculation. The remaining discrepancy then boils down to the relation between the B-meson shape function and the soft fragmentation functioñ s D . It would be interesting to directly computes (2) D from its operator definition in the future, in order to figure out its difference with the shape function.

Numerical results
The contribution of the soft function to the differential cross section is through the factorization formula in the soft limit [33,59,60] dσ N dM d cos θ ∝ ij=qq,gg Tr H ij ln M 2 µ 2 , β, cos θ, µ s ij ln whereN = Ne γ E with N = M/Λ being the moment variable, and H ij are the hard functions, whose matrix elements at the leading order are given by In order to assess the numerical impact of the NNLO correction to the soft function, we define the following quantities: with n = 0, 1, 2 denoting the LO, NLO and NNLO soft contributions, respectively. Here, µ def is the default soft scale, and we exploit two kind of choices for it: µ def,1 = Λ (which corresponds to the choice made in [31]) and µ def,2 = Λ 1 − β 2 cos 2 θ (which was found to be a better choice in [32]). We then define the ratios to quantify the relative size of the NLO and NNLO soft corrections with respect to the LO contribution. In reality, the strong coupling α s should also be evaluated at the soft scale µ, and hence depends on the moment variable N. Here for illustration purposes, we fix its value at α s = 0.118. This does not affect the qualitative behaviors of the soft corrections shown below.
In Figure 8, we show the relative soft corrections as a function of β for µ = µ def , with the two choices of µ def . It can be seen that in the qq channel, the soft function is not very sensitive to the choice of the default soft scale, and the NNLO correction stays below 5% in the whole range of β. On the other hand, the behavior of the gg channel soft function is rather different. With the choice µ def = Λ, both the NLO and NNLO corrections become very large in the boosted limit β → 1, where the NNLO correction can reach about 28%. When changing to the choice µ def = Λ 1 − β 2 cos 2 θ, the corrections are well under control in the boosted region, where the NNLO correction is only about 8%. These findings are in coincidence with the discussions in [32], where the second choice was identified as the better choice, especially for boosted top quark pair production.
In Figure 9, we show the relative soft corrections as a function of µ/µ def for 4 values of β: β = 0.1 (the threshold region), β = 0.5 (the dominant region for the total cross section), β = 0.9 (the boosted region) and β = 0.99 (the ultra-boosted region). For β = 0.1, the two choices of the default soft scale are almost indistinguishable as the solid curves and the dashed curves overlap with each other in the first row of Figure 9. When β gradually becomes larger, the two choices start to differ. We observe that the second option is a good choice in the whole range of β, in the sense that the soft corrections remain small when the soft scale is varied around the default value. This is particularly true for the gg channel, where the t-and u-channel propagators in the tree-level hard function push the average value of cos 2 θ towards unity when the top quarks are highly boosted. In that case the effective soft scale is much smaller than Λ, and is better modeled by the θ-dependent function µ def,2 = Λ 1 − β 2 cos 2 θ. While the above findings have been advocated in [32] from studying the massless soft function, our analysis using the massive soft function provides more comprehensive information on the behavior of the soft corrections in the whole range of phase space.

Conclusion and outlook
In this paper we have calculated the threshold soft function for top quark pair production at the NNLO. We used integration-by-parts identities to reduce the double-real contributions to 23 master integrals, and the virtual-real contributions to 36 master integrals. We then employed the method of differential equations to solve for the master integrals to arbitrary orders in the dimensional regulator ǫ. Our final results are fully analytic, and can be entirely written in terms of GPLs, which makes it efficient for numerical evaluation. Our result represents the first ever NNLO soft function for processes involving a non-trivial color structure and two massive partons with full velocity dependence. Due to the complicated color structure, the renormalized soft function calculated in this paper is a 2×2 (3×3) matrix for qq → tt (gg → tt) in color space. The scale-dependence of the soft function is in full agreement with the two-loop anomalous dimension matrix calculated in [34,35,48] from virtual amplitudes, as expected from the cancellation of infrared singularities. However, the set up of the calculation for the soft function is different from the calculation for virtual amplitudes, and therefore represents an independent confirmation of previous results on the two-loop anomalous dimensions. We find that the previously calculated three-parton correlation function f 2 in the anomalous dimension comes entirely from virtual-real diagrams involving the two massive Wilson lines. This suggests a Coulomb/Glauber origin of the threeparton correlation. It would be interesting to investigate about this in the future.
Our result contains full velocity dependence of the massive partons, which generalizes previous results in more restricted kinematic configurations to fully generic configurations. We have checked that in the limit β → 0, our result reproduces the corresponding results for color singlet/octet production. In the high energy limit β → 1, our result exhibits the expected factorization property of mass logarithms, which leads to a consistent extraction of the soft fragmentation function. We also find full agreement with the NNLO massless soft function in [50], up to the three-parton virtual-real contributions not calculated in that paper.
Our result is an important ingredient in the resummation of threshold logarithms for top quark pair production beyond the NNLL accuracy. It is interesting to study its phenomenological implications for the LHC experiments and future high energy colliders, following the framework of [31,32]. For this purpose, we have studied the numerical impact of the NNLO corrections to the soft function. We find that using a θ-dependent default choice for the soft scale µ def = Λ 1 − β 2 cos 2 θ makes the perturbative behavior well under-control in the whole phase space, from production at rest to highly-boosted production. This is in accordance with the findings of [32]. To achieve the next-to-next-to-next-to-leading logarithmic accuracy for resummation, one also needs the two-loop hard function and the three-loop soft anomalous dimension matrix. The two-loop hard function can be extracted from the virtual amplitudes calculated in [36]. Concerning the three-loop soft anomalous dimensions for multi-leg scattering processes, the result in the purely massless case has been obtained in [45], but the result in the massive case is still lacking and deserves future investigations.
Finally, while our formalism is generic, in practical calculation we have used the fact that the top quark and the anti-top quark are back-to-back in the partonic center-of-mass frame. It would be interesting to consider the more general case, where the tt pair is recoiled by additional particles. This is relevant to the production processes of top quark pairs associated with a Higgs boson or an electroweak gauge boson. At the NLO, the general result of [33] can be utilized, and has been successfully applied to ttH production in [63,64], to ttW production in [91,92], and to ttZ production in [93]. It remains open whether a similar general result can be derived at the NNLO, which we leave for future work.