NNLO soft function for top quark pair production at small transverse momentum

The cross section for top quark pair production factorizes at small transverse momentum of the heavy quark pair, qT. One of the key ingredients that appears in the factorization formula is the soft function, which mediates soft gluon exchanges between particles and gives rise to colour correlations. We present the complete result for the small-qT soft function at the next-to-next-to-leading order. This is the last missing element needed to calculate the NNLO cross section for top quark pair production by means of the qT slicing method. In order to evaluate divergent integrals appearing in the calculation, we develop methods based on sector decomposition and differential equations. We present an extensive validation of our framework. In particular, we recover results predicted by the renormalization group, which constitutes a direct demonstration of validity of the small-qT factorization at NNLO. We provide complete results for the real and imaginary part of the soft function, which are ready for application in the calculation of the tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t} $$\end{document} cross section at NNLO.


JHEP10(2018)201 1 Introduction
One of the most important classes of measurements studied at the Large Hadron Collider (LHC) concern processes which involve the production of the top quark. Of particular interest is top-anti-top production, which is relevant both in studies of properties of the Standard Model, as well as in searches for new physics, as it forms significant backgrounds to many signatures [1].
In the Standard Model, the top quark is the heaviest particle and it does not form bound states, but decays immediately to the W boson and the bottom quark. As the heaviest quark, it also couples most strongly to the Higgs boson, and therefore plays an important role in the mechanism of electroweak symmetry breaking. Top quark pair production enters into pole mass extraction as well as determination of gluon PDFs. Precise predictions for top production allow one also to study rare decays, like those happening through flavourchanging neutral currents, which are predicted to be very small in the Standard Model. The value of the top mass also has an impact on the question of stability of the vacuum [2].
The experiments at the LHC have observed millions of top quarks and more will be detected in the upcoming runs, including the one with high luminosity. The increasing accuracy of experimental data from the LHC is in many cases superior to that of theoretical predictions and this presents a challenge for theory.
The cross section for top pair production is currently known up to next-to-next-toleading order (NNLO) in Quantum Chromodynamics (QCD), both total and differential [3][4][5][6][7]. This, single, complete, NNLO result has been obtained with the approach based on STRIPPER [8][9][10]. There exist also several partial results including NNLO corrections to the off-diagonal channels obtained with the q T slicing method [11], leading-colour NNLO correction to the qq channel calculated with the antenna subtraction [12], as well as approximate NNLO correction including semi-leptonic decays in the narrow-width approximation [13]. Besides, the electroweak corrections for top pair production are known up to NLO [14,15], and they were also combined with the NNLO QCD corrections [16].
In addition to the fixed-order results, a number of resummed cross sections for our process of interest have been calculated. In particular, soft gluons have been resummed at threshold up to next-to-next-to-leading logarithmic (NNLL) accuracy [17][18][19][20][21]. Also, combined resummation of soft and small-mass logarithms at NNLL has been performed [22,23]. Soft and Coulomb gluons have also been resummed simultaneously [24]. Finally, the top quark production cross section has been resummed in the small-q T limit up to NNLL [25][26][27].
Given the complexity of the NNLO calculation for top pair production, a second, independent result is highly desirable. One of the most promising methods that could be used to obtain it is the q T slicing approach [28].
Consider the process h 1 + h 2 → F (q T ) + X, where two hadrons, h 1 and h 2 , collide and produce an object F , which is registered in a detector, together with an undetected QCD radiation X. Then, the cross section at order N m LO can be written as a sum of two JHEP10(2018)201 each of which is separately finite. The advantage of this approach is that the second term in the above equation, which represents resolved emissions, is required only at the N m−1 LO accuracy, and it is already known in most relevant cases. On the contrary, the first term in eq. (1.1), which combines virtual and unresolved real corrections, is usually unknown. However, it is needed only in the small-q T approximation.
In order to calculate the latter, we can use the Soft Collinear Effective Theory (SCET) [29], in which the cross section factorizes at small q T . SCET is an effective theory derived from QCD by expanding diagrams around low-energy scales related to soft and collinear particle emissions. It is based on the strategy of regions [30] and it leads to representing a single QCD field by a set of fields corresponding to collinear, anti-collinear and soft radiation. The hard degrees of freedom are integrated out into Wilson coefficients, which are then used to adjust couplings of the effective theory. The new fields decouple in the Lagrangian and this separation largely facilitates proofs of factorization theorems.
One of such factorizations [25] lies at the basis of the formalism used in our calculation, and the only missing piece needed to use it to evaluate the cross section for top pair production at the next-to-next-to-leading order is the NNLO, small-q T soft function. The result for the latter is presented in this work.
Our result, together with the framework and tools developed to obtain it, form key elements of an alternative calculation of the complete NNLO cross section for the top pair production. They also make up an essential step towards extending the q T slicing method to N 3 LO.
We note that our calculation shares many features with that of the NNLO soft function for top pair production in the threshold limit [31,32]. However, the result for the latter is not of direct use in the context of q T slicing.
The paper is organized as follows. In section 2, we introduce, in detail, all concepts relevant for our calculation. In particular, we define the variables, discuss the factorization and define the small-q T soft function. We elaborate on the colour algebra and introduce the idea of multiplicative renormalization of the soft function. In section 3, we present the LO and NLO soft function, the latter up to the order , which is necessary for renormalization of the NNLO soft function. Section 4 is devoted to a detailed description of the calculation of the bare NNLO soft function. We discuss diagrams which have to be included and we explain the methods which we developed to evaluate all the divergent integrals. Finally, in section 5, we present the results for the complete, small-q T NNLO soft function up to the order 0 . There, we also validate our framework by comparing the results from direct calculation to predictions from the renormalization group, and by comparing the results for a sub-class of NNLO graphs obtained with two different methods. In that section, we also show the results for the NNLO soft function after renormalization. Our findings are summarized in section 6.

4)
where q T is the transverse momentum of the tt pair, y is its rapidity, and m t is the top quark mass. The small-transverse momentum limit is defined aŝ We carry out our calculations consistently in d = 4 − 2 dimensions. It is convenient to introduce the following vectors ≡ k −n µ 2 + k + n µ 2 + k µ ⊥ . (2.8) Notice that k ⊥ is a d-vector, although with pure d − 2-dimensional transverse part. The momenta of the incoming, p 1 and p 2 , and the outgoing, p 3 and p 4 , partons, can be written as where k µ i is a residual momentum which scales like the soft mode k µ i ∼ λ = q T /M 1. The total d-momentum of the tt pair reads q = p 3 + p 4 . (2.10)

JHEP10(2018)201
We also introduce The variable β t is related to the relative velocity of the top and anti-top (or, equivalently, velocity of the top quark in the tt rest frame) [33] | v t − vt| = 2β t . (2.12) We will also work in coordinate space, where the position of the tt pair is given by where x ⊥ is a d-dimensional vector with purely transverse part and the length x 2 T = −x 2 ⊥ . It is useful to express its norm through the following logarithm (2.14) The soft function is invariant with respect to rescalings v i → κ i v i , where κ i s are arbitrary constants. It turns out to be convenient to use this property and redefine the vectors v 3 and v 4 with slightly different normalizatioñ Unless stated otherwise, we shall use the following parametrizations for the d-vectors , sin θ 1 sin θ 2 , sin θ 1 cos θ 2 , cos θ 1 ) , where k and l are the momenta of the two soft gluons radiated at NNLO. The scalar products involving the above vectors read We see, in particular, that At small transverse momenta of the top quark pair, q T , the cross section factorizes according to the formula [25] which is a convolution of the beam functions, B q/h i , B µρ g/h i , the hard functions, H qq , H µνρσ gg , the soft functions, W qq , W gg , and the anomaly exponents F qq (x 2 T , µ), F gg (x 2 T , µ). Above, y corresponds to the rapidity of the top-anti-top quark system, defined in eq. (2.4), θ to the scattering angle of the top quark in the tt rest frame, and φ is the relative azimuthal angle between x ⊥ an v 3 .
We note that the gluon beam functions depend, in general, on the vector x ⊥ , which means that they are sensitive both to its length and to its azimuthal position in the transverse plane. On the contrary, the quark beam functions depend only on the magnitude of the transverse position vector.
The functions appearing in eq. (2.30) capture contributions of gluon emissions from different regions of phase space. In the light-cone parametrization, k µ = (k + , k − , k ⊥ ), the momenta corresponding to each function scale as (2.31) The beam functions are process-independent and they are currently known up to NNLO [34,35]. The hard and the soft functions are not universal, hence, they have

JHEP10(2018)201
to be calculated on a process-by-process basis. The hard function can be extracted from refs. [36,37]. The small-q T soft function has been calculated up to NLO in refs. [25,26].
Each function defined on the right hand side of eq. (2.30) is separately divergent when calculated directly from diagrammatic definitions, like the ones that shall be discussed in section 4. These divergencies correspond to the soft and collinear limits and they must cancel between the hard, soft and beam functions, as the entire cross section has to be finite.
It turns out to be useful to remove divergences also at the level of the functions entering the factorization formula (2.30). This can be achieved by the procedure of multiplicative renormalization, and it will be discussed in detail, for the case of the soft function, in section 2.8. As usual in the procedure of renormalization, the renormalized object acquires dependence on an arbitrary parameter µ, which has to vanish at the level of the cross section. This implies that the renormalized versions of the functions entering the factorization formula (2.30) must satisfy certain evolution equations that govern their µ dependence.

The soft function
The general definition for the position-space, small-q T soft function entering the factorization formula (2.30) reads [25] W whereT and T represent time and anti-time ordering [38]. The normalization factors are (2.33) where N is the number of colours, hence, in QCD, N = 3. It turns out to be convenient to insert the sum over all final states (both discrete and continuous quantum numbers) into the definition (2.32) and obtain where P ⊥ is the total transverse momentum carried by gluons or massless quarks (which recoil against the tt system) in the state |X . The amplitude X|T [O s (0)]|0 is shown in figure 1. The operator O s is defined as a product of Wilson lines, which mediate soft gluon exchanges between particles and give rise to colour correlations. The Wilson line for a particle in representation R, moving along the straight line with d-momentum n i between the points x + an i and x + bn i , is defined as [39]  where P denotes path ordering. The indices a i and b i arise since the Wilson line is an operator in colour space and it accounts for colour evolution of a particle due to gluon emissions. The operator O s for top quark pair production in the qq channel is given by (2.38)

Azimuthal averaging
As we see in eq. (2.30), the factorization formula involves integration over φ, the azimuthal angle between x ⊥ an v 3 . If we restrict ourselves to the NNLO cross section for top quark pair production, the NNLO soft function will enter the factorization formula (2.30) multiplied by the leading order hard function and the leading order beam functions, both of which are Lorentz scalars. Hence, the integration over the azimuthal angle can be pulled directly to the soft function. This motivates us to define the averaged soft function with dΩ d−3 being a rotationally-invariant measure on the unit (d − 3)-sphere, defined recursively as (2.43) Finally, we mention that in order to compute the complete NNLO cross section for tt production in the gg channel, one also needs a version of the NLO soft function which is averaged azimuthally together with the x µ ⊥ x ν ⊥ Lorentz structure which comes from the NLO gluon beam function [35]. This component can be obtained using the unaveraged NLO soft function of ref. [26].

Colour space
The soft function discussed so far is an abstract operator in colour space. It turns out to be useful to represent it as a matrix, with the elements We choose the basis vectors, |c I = c iī I {a} , following ref. [31]. In the qq channel, they read and in the gg channel The inner product is defined as a sum over all colour indices (2.47) The above basis vectors are orthogonal but not orthonormal [31]. The NLO and NNLO contributions to the soft function can be represented as where w are the colour matrices defined as

JHEP10(2018)201
The Hermitian, colour operators T i satisfy the following relations [40] T where 1 is the identity operator, a i ∈ {q,q, g}, and C g = C A , C q = Cq = C F . The operators T i act on vectors in the colour space as follows The matrix elements of the i th parton operator, T i , are given by (T c i ) c 1 c 2 = if c 1 cc 2 , for an initial-state and final-state gluon, , for an initial-state quark (anti-quark).

Fourier transform
In practice, it is easier to calculate the soft function in momentum space. The relevant scalar integrals, which appear in eqs. (2.48) and (2.49), can be transferred to momentum space by means of the Fourier transform where {i} represents the two, three or four particle indices. The structure of the I {i} (x ⊥ ) integrals, for all types of diagrams encountered at NNLO, is where the ellipsis depend on details of the graph, i.e. the direction of the Wilson lines and the soft parton emissions, and p denotes the total momentum of the latter. By plugging eq. (2.54) to eq. (2.53a), we obtaiñ And by applying the azimuthal averaging of eq. (2.39) we get In the above, we used rotational invariance of the measures and the fact that the order of the integrations can be changed. This leads to a reinterpretation of the angles in dΩ d−3 as

JHEP10(2018)201
the azimuthal angles of q ⊥ and that allows one to go from (2.56) to (2.57) with the help of the identity (2.58) Now, we rescale the momenta associated to the emissions by q T , in particularp µ = q T p µ , and get Above, the overall power of q T has been denoted as r. It is a sum of the contribution from the d − 2-dimensional delta function and a genuine contribution from the graph part. The rescaling factorizes all the q T dependence. As a consequence, all integrals in momentum space are proportional to the factor 1/q r T , with the power r, which depends on the order of perturbative expansion. Hence, the transformation from momentum to position space, by means of the Fourier Transform (FT), will amount to multiplication by a factor given by the compact expression (2.60) At NLO, r = 2 + α, and at NNLO, r = 2 + 2α + 2 for double-cut diagrams and r = 2 + α + 2 , for single-cut diagrams, where α is the analytic regulator discussed in the next section. An important feature of (2.60), both at NLO and NNLO, is that its expansion begins at order 1/ Hence, to obtain a result in position space at the order 0 , one needs to calculate the momentum-space soft function up to 1 .

Analytic regulator
The phase space integralsĨ {i} turn out to be divergent not only when the gluons become soft, but also in the limit where the light-cone components of gluons' momenta tend to zero or infinity. Through the relation this limit occurs when the gluon rapidity y g → ±∞. That is why, these are called "rapidity divergencies". They arise because, in SCET, we approximate the full QCD Feynman integrals following the expansion by regions. Yet, we integrate each expression over the full phase space of gluons' momenta. We note that this does not give rise to double counting [29]. Nevertheless, it forces us to introduce another regulator to handle the integrals. The reason why contributions from different regions do not overlap even when integrated from −∞ to +∞ is that integrals in each region depend on a single scale and JHEP10(2018)201 expanding them further (e.g. soft integrals in the collinear region) leads to scaleless integrals, which vanish.
In our calculation, we chose to adopt the analytic regulator prescription of ref. [41], which amounts to the following replacement of the integration measure where ν is a free parameter introduced for dimensional reasons (an analogue of µ in dimensional regularization) and we denote δ + (k 2 ) = δ(k 2 )θ(k 0 ). The regulator α becomes necessary at intermediate stages of the calculation. Since rapidity divergencies do not appear in full QCD, the result for the complete cross section is finite in the limit α → 0. In the case of the analytic regulator, the poles in α cancel not only at the level of the cross section but even at the level of the soft function. This comes from the fact that the soft function for the Drell-Yan process is equal to one, which implies that the product of the beam functions is α-independent (but not the beam functions themselves, see ref. [35]). As beam functions are universal, the same, α-independent product of the beam functions is used in the process of tt production. This means that the only dependence on α can occur in the soft function. Therefore, all α poles have to vanish within the latter. We shall use this feature as one of validations of our calculation.

Renormalization
The renormalized soft function of the small-q T factorization satisfies the following renormalization group evolution (RGE) equation [42] and γ h iī is defined as a non-Γ cusp part of the full anomalous dimension matrix Γ [19], while γ i is the massless-particle anomalous dimension (and enters RGE equations for beam functions in Drell-Yan and Higgs production [43,44]). To make the notation lighter we shall omit the index iī, keeping in mind that the soft function and the anomalous dimension are different in the qq and gg channels.
The soft anomalous dimension matrix γ s is related to the soft renormalization factor (also a matrix in colour space), Z s , as follows and Z s absorbs all UV divergences such that

JHEP10(2018)201
Each quantity in the above equation has a perturbative expansion, either in the renormalized coupling, a s = α s /(4π), or in the bare coupling, a 0 s = α 0 s /(4π), and the relation between the two is where the MS renormalization constant reads and β 0 is the one-loop coefficient of the QCD β function given in eq. (C.5). Hence, substitution of eq. (2.43) for the bare and renormalized soft function, as well as eqs. (2.68) and (2.69) to eq. (2.64) leads to the following, order-by-order relations The quantities on the l.h.s. are finite in the limit → 0. The Z For notational simplicity, in what follows we will absorb the ξ αs prefactors into the definitions of the bare soft functions and change the notation according to Part (I) on the r.h.s. has only terms singular in , which come from the fact that the Z s factors are defined in the MS scheme. These pole terms have to cancel against the singular terms of part (II), which can be used in the following way: knowing Z where L ⊥ was defined in eq. (2.14). Both the soft function and the anomalous dimension have perturbative expansions All quantities in the above equations are renormalized and they are defined in d dimensions.
Because the anomalous dimension matrix starts at the order a s , eq. (2.75) can be solved iteratively. We just need to plug eqs. (2.76) and (2.77) into eq. (2.75) and remember that the renormalized, 4-dimensional coupling a s also depends on ln µ, which implies where we used the expansion of the QCD β function given in eq. (C.4).
After collecting terms at each order, we arrive at the following differential equations Hence, knowing the soft anomalous dimension to order a 2 s , and the soft function to order a 1 s allows one to determine all pieces of the soft function at order a 2 s except the constant, i.e. L ⊥ -independent term. Specifically, at order a 2 s , we get

NLO soft function
The leading order soft function corresponds to the case without radiation. It is given by the following, constant matrices [25,31]   respectively for the qq and gg channels. Because the LO soft function is not divergent, the above result corresponds both to the bare and the renormalized case. Hence, we see that Z The next-to-leading order, bare soft function, is given by where the colour structure is encoded in the w (1) ij matrices defined as comes from renormalization of the strong coupling, see eq. (2.68), and the remaining prefactors arise from the Fourier transform and azimuthal averaging. Only the real-type diagrams contribute to the NLO soft function as the virtual graphs are scaleless and vanish.
The integration measure can be written as d d k = dk + dk − d d−2 k ⊥ and the k + and k − components are integrated from minus to plus infinity. However, the phase space of integration of these light-cone momenta is restricted by δ(k 2 ) θ(k 0 ) appearing in eq. (3.4). This can be easily seen by rewriting the above condition as The delta function fixes k + k − = k 2 T > 0, hence the light-cone components must be both positive or negative. And the theta function, chooses them to be both positive. All in all, the integration over dk + dk − is restricted to the line depicted in figure 3 (left).

JHEP10(2018)201
The NLO soft function has been calculated up to order 0 in refs. [25,26]. However, as can be seen from eq. (2.74), to renormalize the NNLO soft function, we need to know the NLO soft function up to the order 1 . This order has been calculated in ref. [45], yet, using a slightly different definition of azimuthal averaging. We have adjusted this result to our definition, as well as fully cross checked it using a sector decomposition-based approach (described in detail in section 4.4), finding a perfect agreement.
Hence, the final result for the position-space NLO soft function up to order 1 reads  where and the result for J 34 can be expressed in terms of ordinary polylogarithms.

NNLO soft function: methods of calculation
To calculate the next-to-next-to leading order contribution to the bare soft function, one needs to sum several groups of diagrams, each multiplied by a proper colour factor. The relevant master formula can be derived directly using definitions given in eqs. We shall now discuss groups of diagrams contributing to each of the terms in the above equation. Then, we will describe the methods used to calculate them.
Double-cut diagrams. The gluons in double-cut diagrams can connect two, three or four Wilson lines. Amongst the two-Wilson-line graphs, one can distinguish a special class of bubble diagrams, depicted in figure 4. Besides gluons, they may also involve quarks and, as we work in the Feynman gauge, also the ghost bubble. Even though they belong to the class of two-cut diagrams, they are easier to calculate than most of the graphs in this group. The integrals corresponding to the bubble diagrams, and methods applied to evaluate them, shall be discussed in detail in section 4.2. The bubble graphs come with colour matrices which are identical to those appearing in the NLO soft function, i.e. w (1) , given in appendix A. The second group consists of non-bubble graphs in which the gluons attach to only two distinct Wilson lines: i an j. These graphs are shown in figure 5. We see that they include both abelian and non-abelian structures.
The above observations were also made in the related calculation of the NNLO soft function for top pair production in the threshold limit [31,32].
The abelian graphs depicted in figure 6 constitute the only non-vanishing three-Wilsonline contribution in the group of double-cut diagrams. This is because the non-abelian, three-Wilson-line graphs cancel when summed over colour structures, which can be understood by analyzing the graph depicted in figure 7. Since it is a double-cut diagram, the corresponding integral is a real-valued function. And it is multiplied by the colour factor if abc T a i T b j T c k . The complete soft function receives also a contribution from a diagram which is a complex conjugate of figure 7 and the complex conjugation only affects the colour factor, turning it into −if abc T a i T b j T c k . Hence, the diagram of figure 7 and its complex conjugate cancel, and the non-abelian graphs with three distinct Wilson lines connected by the gluons do not contribute to the soft function.
We emphasize that this happens only because the integrals are real functions, as they originate from double-cut diagrams. This property will not be true for single-cut, nonabelian diagrams, and we will see that the corresponding contributions do not vanish when summed over all diagrams.   The complete expression for the NNLO soft function derived using the diagrammatic approach described above can be obtained alternatively by taking the soft limit of the relevant matrix elements in QCD.
In the case of the gg final state, in the limit k, l → 0, the squared matrix element factorizes as [10] which we had used to construct the NLO soft function integrals in eq. (3.4), and p i are the d-momenta of the external partons. We note that the functions S ij (k) and S ij (k, l) are invariant with respect to rescalings of the momenta of external particles of the Born process. Therefore, they can be expressed in terms of velocities The function S ij (k, l) can be split into two parts where m i and m j are the masses of the external particles. The first term in eq. (4.5) has been given in ref. [46] and reads In the above equation, the first line comes solely from the gluon and ghost bubble diagrams of figure 4. The second line originates from the C A part of diagrams which do not involve the triple-gluon vertex, that is D 4 and D 5 in figure 5. Finally, the last two lines receive contributions from the non-abelian diagrams D 6 , D 7 and the gauge bubble. The second contribution in eq. (4.5) was derived in ref. [9] and represents additional terms generated by non-vanishing masses. The relevant function is Here, the first and the third term come from the C A part of diagrams which do not involve the triple-gluon vertex, whereas the second term arises due to non-abelian contributions, including the gauge bubble. In the case of the final-state qq-pair, in the limit k, l → 0, the matrix element factorizes as [10] where the function I ij (k, l) has the form The above expressions can be used directly to calculate our soft function of interest by applying the following formula where the case with f = q corresponds to the qq final state while the case with f = g corresponds to the gg final state. We have checked through direct calculation that the expression for the NNLO soft function derived using definitions given in eqs. (2.32)-(2.38), i.e. expanding the Wilson lines and truncating the series at O α s 2 , matches exactly the formula (4.10).
To complete the definition of eq. (4.10), we need to specify the colour matrices w Single-cut diagrams. The single-cut diagrams required for the calculation of the NNLO soft function for top pair production are shown in figure 9. In this class of graphs, one gluon is real and the other runs in a loop. Because of the latter, the integrals from the single-cut graphs produce also an imaginary contribution.

JHEP10(2018)201
As in the case of the double-cuts , also the single-cut term derived using the diagrammatic approach described above can be obtained alternatively by taking the soft limit of the relevant matrix elements in QCD. This results in the formula [47] where while R ij and I ij correspond to the real and imaginary parts from integration over the loop momentum and they were obtained in refs. [47,48]. We notice that a new, antisymmetric colour structure appears in the imaginary part: ijk is defined in eq. (2.50c) and its explicit expressions are given in appendix A.
We have checked through direct calculation, based on definitions given in eqs. (2.32)-(2.38), that we reproduce the single-cut expression of eq. (4.11). With the latter, we can then determine the real-virtual part of our soft function Zero-cut diagrams. Purely virtual, two-loop diagrams, do not involve a cut gluon. Therefore, the measurement function δ (d−2) (f (k, l) − q ⊥ ) does not appear in the corresponding integrals. As a consequence, these integrals are scaleless and vanish in dimensional regularization. Hence, the NNLO soft function for top pair production does not receive contributions from two-loop diagrams.

Symmetries between integrals
Our double-cut, soft function integrals have the general structurẽ where p i are the momenta of external particles. The above integral can depend only on the scalar productions p i ·p j , which are invariant with respect to Lorentz transformations p µ → Λ µ ν p ν . To balance the transformation of the external momenta in the scalar products p i · k and p i · l one needs to transform the gluon momenta k and l with Λ −1 µ ν . This, however, does not leave the integral unchanged for a general Lorentz transformation because of the transverse delta function appearing JHEP10(2018)201 in eq. (4.14). Therefore, the integralĨ({p i , p j }) is invariant only under a subgroup of the Lorentz group which involves: • rescaling of the light-cone components of the momenta k and l compensated by inverse rescaling of the light-cone components of the external momenta, • rotations of the above momenta in the transverse plane such that |k ⊥ + l ⊥ | = q T .
Let us denote Given the above reduced Lorentz symmetry, and the fact that the result can only depend on external momenta, we conclude, that our integral must be a function of p i * p j and p i,⊥ · p j,⊥ or, equivalently, p i * p j and p i · p j .
As is clear from eqs. (4.3)-(4.9), the integrand is invariant under the rescaling of p 3 = λ 3 p 3 and p 4 = λ 4 p 4 , for arbitrary λ, which means that the result can only be a function of ratios of the above scalar products.
Integrals which do not exhibit α poles are, in addition, invariant under rescaling of p 1 . But, here, one cannot form a variable of the type p 1 · p 3 / p 2 1 p 2 3 because p 1 is massless. As a consequence, the only possibility is p 3 * p 3 /p 2 3 = p 4 * p 4 /p 2 4 . This is why some abelian integrals are equal. In the cases of the integrals with α poles, rescaling of p 1 is broken by the analytic regulator: I(λp 1 ) = λ −α I(p 1 ).

Differential equations approach and bubble diagrams
With an exception of the propagators introduced to regularize rapidity divergences, the quark, gluon and ghost bubble, as well as part of the triple gluon vertex diagrams depicted in figure 5, which we shall call "tadpole", depend only on the momenta k and k + l. This feature allows one to first integrate over k + l and then solve the integral over k with help of the differential equation approach [49,50]. The integration can be performed mostly analytically, with an exception of a few one-dimensional integrals which we integrate by numerical methods.

Vacuum polarization tensor
We consider the process depicted in figure 10 (a), where the particle running in the loop can be a quark, a gluon or a ghost. The corresponding vacuum polarization integrals read where C f = C A for the gluon bubble and the tadpole, while C f = −T F for the quark bubble. The numerator N µν (k, p) depends on the particle in the loop. For the quark loop we have for the gluon+ghost loop and for the tadpole The directions of momenta are indicated with arrows in figure 10 (a). All propagators are assumed to be defined with the +i prescription. The most generic self-energy tensor that can be formed out the d-vectors p and n, and the metric tensor g µν is ImΠ µν (α) = T 00 g µν + T pp p µ p ν + T nn n µ n ν + T pn (n µ p ν + p µ n ν ) . (4.21) The coefficients T ij can be expressed in terms of two scalar integrals, A 0 and B 0 , through the procedure of Passarino-Veltman reduction [51] A (a) (4.23) The integral A 0 can be shown to be scaleless and, therefore, it vanishes. Hence, it turns out that the only two-point integral that we need to determine ImΠ µν is B 0 , and this integral can be calculated exactly by means of Schwinger parametrization. The result takes the form Unitarity allows us to obtain also a version of this function which corresponds to the cut diagram of figure 10

JHEP10(2018)201
Finally, the coefficients in eq. (4.21) take the following forms. For the quark bubble For the gauge (i.e. gluon + ghost) bubble And for the tadpole The function B 0 (α, α) in eqs. (4.27), (4.29) and (4.30) can be given either by eq. (4.24) or by eq. (4.25). In our soft function calculation, we will use the latter, cut version of this two-point integral.
One can check that the results given in this section recover standard expressions for the gluon vacuum polarization tensor in the limit α → 0, where, in particular, the coefficients T nn and T pn vanish.

Soft function integrals
The above results for the quark and gauge bubble can now be embedded in the two-Wilsonline soft function graphs, as depicted in figure 4. The tadpole integrals correspond to the diagram D 6 of figure 5 with the eikonal propagator connecting the two gluons replaced by a pinch. The integrals that appear have the following structure We note that, as in the case of the NLO soft function discussed in section 3, the integration measure can be written as d d k = dk + dk − d d−2 k ⊥ and the k + and k − components are integrated from minus to plus infinity. However, the phase space of the integration over the light-cone momenta is restricted, this time by θ(k 2 ) θ(k 0 ) = θ(k + k − − k 2 T ) θ(k + + k − ). Because of the first theta function, the integration is fixed not to a line, as in the case of the NLO soft function depicted in figure 3 (left), but to the region defined by the condition JHEP10(2018)201 T , which means that k + and k − have to be both positive or negative. The second theta function chooses k + and k − to be both positive. Hence, the d-momentum of the off-shell gluon which appears in the bubble and the tadpole diagrams is integrated over the region depicted in figure 3 (right).
We would like to solve the class of integrals from eq. (4.31) by means of the method of differential equations [49,50], with help of reverse unitarity [52], which allows one to turn delta functions into propagators. What prevents us from direct use of the latter is the θ(k 2 ) function. However, we can trade this theta function for the Dirac delta function at the cost of introducing an extra integration over a spurious mass, m 2 . Namely, we multiply eq. (4.31) by which leads to (4.33) Hence, we obtain Now we can use reverse unitarity to turn delta functions inĪ(m 2 ) into propagators, which leads to the following topologȳ I(a 1 , a 2 , a 3 , a 4 , a 5 , a 6 ) = (4.35) We observe that all poles, hence those of the soft origin, are generated by performing the integral over m 2 in eq. (4.34) while the functionĪ(m 2 ) is finite in the limit → 0.
A set of identities can be derived for the class of integrals defined in eq. (4.35). First of all,Ī(a 1 , a 2 , . . . , a 6 ) obeys the standard integration by parts (IBP) identities with q µ = n µ ,n µ , v µ 3 , v µ 4 , k µ . The set of IBPs consists of five relations. In addition, the topologyĪ exhibits certain redundancies which result in additional identities. One of them comes from the q T -delta propagator
(4.38) The last identity arises from momentum conservation following the discussion of section 2.1 n +n =ṽ 3 +ṽ 4 . (4.39) Multiplying the above by k leads tō I(a 1 , a 2 , a 3 − 1, a 4 , a 5 , a 6 ) +Ī(a 1 , a 2 , a 3 , a 4 − 1, a 5 , a 6 )− 1, a 2 , a 3 , a 4 , a 5 , a 6 ) −Ī(a 1 , a 2 − 1, a 3 , a 4 , a 5 , a 6 ) = 0 . Hence, altogether we have seven identities which we use to reduce all the integrals appearing in the problem to a set of master integrals. While solving the bubble graphs with the above method, it is important to note that the d-vector k µ in eq. (4.35) is now massive. Hence, the parametrization from eq. (2.16) does not hold and it must be replaced with k = (k 0 , . . . , | k| sin θ 1 sin θ 2 , , | k| sin θ 1 cos θ 2 , | k| cos θ 1 ) , Therefore, the relevant inner products now take the form n · k = k 0 − | k| cos θ 1 , (4.43) n · k = k 0 + | k| cos θ 1 , (4.44) v 3 · k = k 0 − β t | k| sin θ 1 cos θ 2 sin θ − β t | k| cos θ 1 cos θ , (4.45) v 4 · k = k 0 + β t | k| sin θ 1 cos θ 2 sin θ + β t | k| cos θ 1 cos θ . We also get that k 0 = 1 2 (n · k +n · k) = 1 2 (ṽ 3 · k +ṽ 4 · k) . (4.48) Since the values of the powers a i which appear in the definitions of the integrals are governed by the powers in the denominator of eq. (4.17) and the powers of the p 2 , and n · p propagators in the expressions for the coefficients T ij , eqs. (4.26)-(4.30), they are the same the for the quark and gauge bubble, as well as for the tadpole integrals. Therefore, obtaining a solution for one type of the bubble, allows us to use it for the other type.

(4.52)
We start by reducing them with help of IdSolver [53] -a C++ implementation of the Laporta algorithm [54], which depends on FORM [55] and Fermat [56]. As a result, we obtain five master integrals, for which we then derive a set of differential equations with respect to the variable β t . The structure of the set is such that the general solutions for the masters can be obtained iteratively as a series in α and .
In fact, we only had to solve two systems of differential equations: one for three masters from the reduction ofĪ 13 and the other for two masters from the reduction ofĪ 34 . The reduction ofĪ 33 leads to master integrals identical to those found from reduction ofĪ 13 . The same masters can be used to determineĪ 23 as discussed above.
The expressions for the boundary integrals, corresponding to β t = 0, are easily calculated through direct integration and this allows us to determine the special solutions. As a last step, we integrate the expressions for the functions from eq. (4.52) over the spurious mass m 2 , following eq. (4.34). Except for a few one-dimensional integrals which appear at order 1 in momentum space, most of the result is given in an analytic form.

Real-virtual diagrams
We now turn to the class of single-cut diagrams. As shown in figure 9, these diagrams involve a gluon loop as well as a real gluon. Following ref. [47], we introduce the tree-level, UV-renormalized, one-loop soft currents where e µ i is the eikonal propagator defined in eq. (4.12). The function g ij , which is symmetric under the exchange i ↔ j, has been obtained in a concise form in ref. [48]. The soft current J µ a corresponds to the sum of all parts of diagrams shown in figure 9 to the left of the cut. It can therefore be used directly to construct the single-cut contribution to our soft function of interest by attaching the gluon with momentum p to the Wilson lines in a Born-level amplitude. As a result, we obtain where [48] g (1) ij ( , p, p i , p j ) = − ij , were derived in refs. [47,48], and they depend solely on the rescaling-invariant variables . (4.57) From eq. (4.56), and taking into account the expansion of the Fourier Transform coefficient of eq. (2.61), we observe that the single-cut contributions to the position-space NNLO soft function start at the order 1/ 3 for the real part and 1/ 2 for the imaginary part. As we shall see in section 5, the 1/ 3 term cancels between the single and the double-cut contributions.

Diagrammatic configurations
A range of different configurations of diagrams contribute to eq. been calculated in ref. [57]. Connecting this current to a massive Wilson line leads to an expression which can be integrated analytically.
Another interesting subclass of single-cut diagrams is formed by two-Wilson line configurations with gluons attached to two massless and one massive leg (e.g. 131). These integrals correspond to case 1 of ref. [47] and they exhibit rapidity singularities. Finally, we also need to include diagrams with three distinct Wilson lines, which corresponds to case 3 of ref. [47].
Altogether, the single-cut part of the soft function receives contributions from ten independent two-Wilson line and ten independent three-Wilson line integrals. All the other integrals can be obtained through symmetry relations.
In particular, for the two-Wilson line diagrams we observe the symmetry D i,ji = D i,ij , where the comma corresponds to the cut in the diagram. Swapping i ↔ j on one side of the cut has no effect as the colour and the kinematic parts both produce minus signs which balance each other. One can also show that D ij,i = D * i,ij , which has an important consequence as it implies that the two-Wilson line, single-cut diagrams are purely real.
We also find symmetries for the three-Wilson line diagrams. First of all, swapping particles on one side of the cut has similar effect to the one described above for the two-Wilson-line case, namely D k,ji = D k,ij , which, again, arises because changes in signs in the colour and the kinematic part balance each other. The most important relation, however, reads D ij,k = −D * k,ij , and it means that the three-Wilson line, single-cut diagrams are purely imaginary. In fact, this class of diagrams, constitutes the only source of the imaginary part of the entire NNLO soft function for top pair production.

Three-particle diagrams with two massless and one massive Wilson lines
Let us consider a special case of three-Wilson-line diagrams which involve two massless and one massive particle. Two such diagrams are depicted in figure 11. Using the notation introduced above, we can write expressions corresponding to those two diagrams as 12 , (4.58) Hence, their sum reads 12 − g (1) * 12 . (4.60)

JHEP10(2018)201
We see that the result is purely imaginary and that, contrary to the case of double-cut diagrams, it does not vanish as, here, the antisymmetry of the colour factor under the exchange 1 ↔ 2 is compensated by the antisymmetry of the kinematic part.
We note that our case is different from the one of ref. [58], which considers infrared singularities of QCD amplitudes and where it is argued that all three-particle structures with two massless and one massive Wilson lines must vanish. The latter happens because ref. [58] discusses amplitudes, i.e. objects in which all Wilson lines originate from the same vertex and all soft gluons are virtual. On the contrary, in our real-virtual diagrams, the Wilson lines meet at two vertices and one gluon is cut.
To understand better why the contribution from the diagrams of figure 11 does not vanish, it is useful to perform an analysis similar to the one of ref. [59]. For that purpose, we use the soft current for massless particles with momenta p 1,2 = (1, 0, 0, ±1), which reads [57] J a(1) Embedding the above in the integral over p gives Let us now apply the change of variables, motivated by ref. [59] Our integral becomes Hence, we see that the above integral would exhibit scaling w.r.t. ξ, and vanish, if only there was no α regulator. However, without the regulator, the integral is divergent. Hence, we conclude that the contribution from tree-particle graphs with two massless and one massive Wilson lines does not vanish. This result only affects the imaginary part of the soft function. The tree-particle graphs of the type of 123 do not contribute to the real part as there will always be a complex-conjugate diagram with opposite sign due to colour operator.
The above conclusion does not invalidate the analysis of ref. [59], where it is claimed that the massless-massless-massive diagrams vanish because of the scaling property (4.63).

JHEP10(2018)201
The key difference between our case and the one discussed there is that, in ref. [59], purely virtual diagrams are considered. In those diagrams, rapidity divergences are regulated by dimensional regularization [41]. Hence, no α regulator is required and the ξ scaling of eq. (4.63) holds.
The key element of our calculation, which prevents the diagram 123 from vanishing, is the transverse delta function. As explained in ref. [41] without this function, integration over the transverse momentum provides a factor k − − which regularizes rapidity divergences, hence there is no need for the regulator α. However, when the transverse delta is present, it fixes q ⊥ to some external value and the integration over q ⊥ does not provide a regulator of light-cone singularities. Hence, we need to introduce the analytic regulator α.

Method of integration
To evaluate the integrals in eq. (4.55) we proceed through the following set of steps. We start from integration over the p T and p − components, which is straightforward as, in the process, we use the two delta functions δ(p 2 ) and δ(p T − 1). Then, we are left with a nontrivial integration over p + and cos φ , where φ is the azimuthal angle between v 3⊥ and p ⊥ . The remaining angular variables do not appear in the integrand and they only produce a surface term.
To deal with the nontrivial 2-dimensional integrals, we remap the variables p + and φ to a unit hypercube The second transformation is introduced for efficiency reasons as it eliminates integrable singularity in the azimuthal integration. The first transformation compactifies the plus component of the gluon momentum. This is useful as our integrals are in general divergent when integrated over p + . The divergence comes from rapidity singularities and that is why we have introduced the analytic regulator in eq. (4.55). To perform the integral over x, we use the Laurent expansion 66) which makes the divergences explicit by turning them in α poles. All the coefficients of α expansion are finite and take forms of one and two-dimensional integrals which we perform numerically with help of the Cuba package [60,61]. The integration is fast, hence, we are able to achieve arbitrary accuracy.

Sector decomposition approach and real-real diagrams
The double-cut integrals take the following form Two issues arise when one attempts to evaluate them. First of all, they are divergent when integrated over a subset of variables, and the pattern of divergencies is complex, with JHEP10(2018)201 many overlapping singularities. Secondly, the remaining part of the integration, where divergencies do not appear, consists of complicated azimuthal integrals and care is needed to perform them efficiently. Because of these two separate challenges, it is convenient to factorize the graphdependent part as Graph(n i , k, l) = (Graph(n i , k, l)| βt=0 ) boundary part Graph(n i , k, l) Graph(n i , k, l)| βt=0 weight ≡ I(n i , k, l) W(n i , k, l) .
(4.68) For any graph, the boundary part is not only independent of β t but also of all the angles except the angle between transverse components of gluons d-momenta θ 1 = (k ⊥ , l ⊥ ). In other words, all divergences are present already in the boundary integrals and the remaining integration over the angular variables which appear in the weight is finite.
Our strategy of evaluating the double-cut integrals will therefore consist of two elements: (i) use of sector decomposition to disentangle overlapping singularities and cast them into a set of α and poles, (ii) supplementing the sector-decomposed integrals with carefully parameterized weights and integrating them numerically with Cuba.

Integration of the on-shell and transverse delta functions
The momenta of the heavy quarks can be written using the following parametrizatioñ wherev 3⊥ is a unit vector in d − 2-dimensional space, cf. eq. (2.16). As discussed in section 4.1, the integrand in eq. (4.67) can only depend on the scalar products between k ⊥ , l ⊥ and v 3⊥ . Due to rotational invariance of these scalar products, one can always change the frame in the transverse space such that Since the momenta of the incoming particles are light-like and v 4⊥ = −v 3⊥ , the soft function integrals can only depend on v 2 3⊥ . Therefore, nothing is sensitive to the position of the versorv 3⊥ . Hence, the integral (4.67) can be written as where the surface of unit-sphere, S d−3 , is given by eq. (2.42), and the differential measures read Because we use the parametrization (4.71), all the angles except θ 1 , φ 1 and φ 2 can be integrated trivially. Let us now replace the angle between the transverse components of the two gluons' momenta, θ 1 , by Then, the whole measure takes the form where we used eqs. (2.40), (4.73) and (4.74). After performing the graph-independent integrations over k − , l − and η, where the last quantity gets fixed to the integral (4.67) turns intõ and the delta functions generate the following relations and The last pair of inequalities defines the integration region in the (k T , l T ) plane, which is depicted in figure 12 (left). The first inequality corresponds to the two (red) solid lines, while the second inequality corresponds to the (blue) dotted line.

Divergencies
The non-integrable divergences of the integral (4.77) come from the graph-dependent part.
Since the weight is finite, the divergences are fully determined by the boundary integrand I(n i , k + , l + , k T , l T ). Hence, each integral can in principle diverge due to integration over the four independent variables: k + , l + , k T and l T , and some of the limits are coupled due to the constraint of eq. (4.79). The singularities of the integralĨ correspond to vanishing of the propagators in I, and this can happen in the following situations.

JHEP10(2018)201
To see that the massive-particle propagators can produce divergencies only in the soft region, let us use the parametrization of eqs. (4.71) and the relation (4.78). Together they give Since β t ≤ 1 − m 2 t /ŝ < 1 and −1 ≤ cos θ ≤ 1, the coefficients in front of the components of gluon's d-momentum in the first two terms in eq. (4.83) can never vanish. The collinear region corresponds to k T → 0 and finite k + , and we see that the propagator stays finite in this limit. Similarly, the propagator does not vanish when k + → 0 or k + → ∞, while k T = 0, which is the region of small/large rapidities of the gluon. Hence, the only way to make the above propagator vanish is to send k + and k T to zero simultaneously, and this corresponds to the soft limit. The above proof can be repeated for the other propagators of the outgoing particles.
Exact gluon propagator. Double-cut diagrams with the triple gluon vertex contain the propagator (k + l) 2 in which the momenta of the gluons are commensurate. After the integration over k − , l − and η is performed, this exact propagator reads Within the domain of integration defined in figure 12 (left), the above propagator vanishes when k T → k + k + + l + and l T → l + k + + l + . (4.85) This limit corresponds to the two gluons with momenta k and l becoming collinear. Note that, due to condition (4.79), only one of the two gluons can have vanishing transverse momentum. One notices that the limit (4.85) is different with respect to the cases discusses earlier. While the divergencies of the eikonal propagators happen only at the edges of the integration domain, i.e. 0 or ∞ (endpoint singularities), the exact propagator vanishes inside the integration region, on the manifold defined by eq. (4.85).

Mapping procedure
We would like to transform the integration region in the (k + , l + , k T , l T ) space into a unit hyper-cube in the space of variables {x i }, and, if necessary, split it such that each of the resulting integrals i dx i f ({x i }) has singularities only when one or more variables x i → 0. We achieve the above in the following four steps: Step 1. The transverse variables (k T , l T ) are transformed to new variables (x T , y T ) with the replacements  Figure 13. Manifold where the double-cut integrals with triple-gluon vertex become divergent. It corresponds to the limit of the two gluons k and l becoming collinear to each other. and the inverse transformation reads This results in change of the integration region from the one shown in figure 12 (left) to that of figure 12 (right). The latter corresponds to The collinear divergences now occur in the limits and the soft divergences correspond to (4.90) The limit of two gluons becoming collinear, eq. (4.85), in the new variables happens when The last divergence occurs on a manifold inside the integration region as depicted in figure 13.
Step 2. In order to transform the manifold singularity (4.91) into an endpoint singularity, we split the integration region in the variable y T precisely on the manifold of figure 13 dy TĪ (k + , l + , x T , y T ) where c(k + , l + ) = k + − l + k + + l + .
(4.93) Now, we use two different parametrizations (4.94) and (4.95) Hence, we obtain (4.96a) With the above changes, the soft singularities happen at in both integrals, I u and I d , and for any value ofȳ T . The case in which a gluon is collinear to the incoming parton corresponds toȳ for both integrals. We note that in the above limit, for I u , it is the gluon with momentum l that becomes collinear, whereas, for I d , the limit (4.98) corresponds to the gluon with momentum k being collinear to the beam. Finally, the limit of two gluons becoming collinear to each other, eq. (4.85), in the new variables corresponds to x T → 0,ȳ T → 0 , (4.99) for both I u and I d .
Step 3. We see that the integrals (4.96) can be divergent at both ends in the variablē y T . In order to move all singularities to the limit x i → 0, we split I u and I d atȳ Then, we apply the following transformations

JHEP10(2018)201
At this point, we have the following singularities: (4.103) Step 4. In the last step, we compress the ranges of the k + and l + integrals by the transformation whose inverse reads Hence, finally, our integral is a sum of four contributions (4.106)

Integrating the weight
As a final step, we integrate the weight over the azimuthal angles of v 3 : φ 1 and φ 2 . In order to map the angular variables into a unit hypercube, we first write the integration measure as Then, we represent the remaining measure as [10] where we used the fact that our integrand is independent of the angles φ 3 , . . . , φ 1−2 . The plus distribution is defined as The integral over the weight now reads Figure 14. Quark bubble contribution to the NNLO soft function in the qq channel. The plots correspond to three independent entries of the 2×2 matrix determined for the fixed values of θ = π/4 and L ⊥ = 0. We see perfect agreement between RG prediction and the direct calculation obtained with two independent methods: differential equations (DE) and sector decomposition (SD).

JHEP10(2018)201
To map this integral to a unit hypercube, we use the following changes of variables in the second and the third line, respectively. From these, one gets where we define the prefactor (4.114)

NNLO soft function: results
The bare NNLO soft function for the top pair production has the following structure Because it encodes single and double-soft limit, it exhibits at most 1/ 2 singularity. However, higher order poles, as well as α poles, appear at intermediate stages of the calculation, but they cancel in the final combination. and L ⊥ = 0. We see perfect agreement between RG prediction and the direct calculation obtained with two independent methods: differential equations (DE) and sector decomposition (SD).
We have checked that, in our calculation, all α poles, including /α, as well as 1/ 4 vanish within each colour structure defined in eq. (2.50). As for the 1/ 3 pole, it comes out with a non-zero coefficient in the single-cut and in the double-cut contributions. The value of the 1/ 3 coefficient is however identical, in those two pieces up to a sign. Hence, in the final combination, this pole does not appear in our result. We have demonstrated that, when all contributions to the bare NNLO soft function are included, following eq. (4.1) the soft function indeed shows at most 1/ 2 singularity.
As discussed in section 2.8, the pole part of the soft function, i.e. the functions S (2,−2) (L ⊥ , β t , θ) and S (2,−1) (L ⊥ , β t , θ) defined in eq. (5.1), as well as the L ⊥ -dependent part of the finite contribution, S (0,0) (L ⊥ , β t , θ), can be completely determined from the renormalization group. The only term that has to be obtained through direct calculation is the L ⊥ -independent part of S (2,0) (L ⊥ , β t , θ). In spite of the above, it is worth calculating all the terms appearing in eq. (5.1) and use the redundant ones for cross checks of our framework against RG prediction.
As much as the comparison of the results from direct calculations to the prediction of the renormalization group is valuable, it is limited by the fact that RG misses the constant JHEP10(2018)201 part of S (2,0) (L ⊥ , β t , θ). In our calculation, we are however able to cross-check this missing component as well since part of our soft function, namely the one corresponding to the bubble and tadpole diagrams, can be determined with two completely different methods: differential equations and sector decomposition. Figures 14 and 15 show the results for the quark bubble part of the soft function, respectively in the qq and the gg channel, obtained from the renormalization group and from the direct calculation with differential equations and with sector decomposition. The quark bubble contribution can be singled out from the RG prediction as it is proportional to n f . The plots correspond to three independent entries of the 2 × 2 matrix and six independent entries of the 3 × 3 matrix, determined for the fixed values of θ = π/4 and L ⊥ = 0, and shown as a function of β t . We see perfect agreement between the two sets of points as well as between the points and the RG prediction. The error bars in the result from sector decomposition come from the uncertainties of numerical integration. We have performed similar comparison of the results from differential equations and sector decomposition for the gauge bubble and for the tadpole, in both channels, each time finding an excellent agreement.
This constitutes a very strong validation of our computational framework. It tests all the elements of the sector decomposition method and tools, discussed in section 4.4, up to the finite order in . Even though the bubble graphs are easier to solve analytically than the rest of the double-cut graphs, they pose a serious challenge to the sector decomposition approach due to non-trivial numerators. Therefore, the comparison shown in figures 14 and 15 makes up a highly-nontrivial test of the entire framework used in the calculation of the NNLO soft function.
We are now in a position to present the complete NNLO soft function for top quark pair production at small q T . The results for the real part of the independent entries of the matrices in the qq and gg channel, at fixed values of θ = π/4 and L ⊥ = 0, are given in For the poles, we also show, as lines, predictions from the renormalization group. We see that the two sets of predictions agree perfectly.
For the finite term, we show two sets of points: the one before and the one after renormalization. They differ by a finite function following the formula (2.74). We notice that the difference is often substantial.
As discussed in section 4.3, due to the real-virtual contributions, the NNLO soft function contains also an imaginary part. As follows from eqs. (A.9) and (A.10), discussed in appendix A, the imaginary part of the soft function matrix has only one independent element. The corresponding results in the qq and gg channel are given in figure 18. As in the previous cases, the predictions were obtained for the values of θ = π/4 and L ⊥ = 0, as functions of β t . Like in the case of the real part, we see excellent agreement between the renormalization group predictions and the direct calculation for the pole part.

Summary
A second, independent calculation of the NNLO correction for top quark pair production in proton-proton collisions would be of great value. One of the methods that could be JHEP10(2018)201 used to obtain it, is the q T slicing approach. In this article, we presented the result for the complete NNLO, small-q T soft function for tt production, which forms a significant step on the way towards this goal. Now, all the ingredients needed to construct the NNLO cross section for top pair production valid at small-q T are available.
To obtain our results, we have constructed a framework based on sector decomposition and differential equations. The framework has been extensively validated. In particular, we have checked that all the α poles, including /α, and all the poles beyond 1/ 2 , vanish in the complete result for the NNLO soft function. We also checked that we recover all pieces predicted by the renormalization group, both the pole terms and the L ⊥ -dependent part of the finite term. The strongest validation of our framework and tools comes from finding a perfect agreement between numerical results from the sector decomposition-based framework and analytic results obtained with the method of differential equations, for the graphs involving gauge, ghost or quark bubble.
The NNLO, small-q T soft function can now be used to obtain full tt cross section at NNLO by means of the q T -slicing method, as well as for small-q T resummation at NNLL'.

JHEP10(2018)201
resources of ACC Cyfronet AGH where our numerical calculations were performed on the Prometheus cluster. S.S. would like to acknowledge the Mainz Institute for Theoretical Physics (MITP) for hospitality and support during the completion of this work.

A Colour matrices
In this appendix, we collect all colour matrices that appear in the formulae for the soft function at NLO and NNLO, in the basis defined in eqs. (2.45) and (2.46).
Let us start by determining the set of independent dot products of the operators T i . For the top pair production in the iī channel, we have The colour conservation allows us to write the following set of equations for i, j ∈ {1, . . . , 4}. By combining (A.1) and (A.3), we arrive at the relations Hence, in general, we have four independent structures: For the qq channel, this number is further reduced to three, as all T i · T i are equal, which leads to the relation T 1 · T 2 = T 3 · T 4 . In the case of the qq channel, the three independent matrices w (1) ij , defined in eq. (2.50a), read

JHEP10(2018)201
while, for the gg channel, we have the following four independent w (1) ij matrices w (1) The w (2S) ij matrices, which enter the abelian part of the double-cut contributions to the soft function (4.2), were defined in eq. (2.50b). In the case of the qq channel, they take the following explicit forms In the gg channel, these colour matrices read

JHEP10(2018)201
As for the antisymmetric configurations, w (2A) ijk , there is only one independent matrix per channel, which we choose to be w All the other antisymmetric matrices can be obtained as combinations of the above and w (1) ij , using colour conservation (A.2), together with the relation

B Cusp angles
The amplitudes are functions of Lorentz invariants [62] s ij ≡ 2σ ij p i · p j + i , To this end, we label massless particles with lowercase i, j, . . . and massive ones with capital indices I, J, . . .. The soft anomalous dimension of eq. (2.66) is a function of cusp angles, β ij , β Ij and β IJ , formed by massless and/or massive Wilson lines. They read [62] β ij = ln −2σ ij p i · p j µ 2 (−p 2 i )(−p 2 j ) where we have introduced which comes from regularization of the IR divergences in the effective theory by taking massless partons slightly off-shell (−p 2 i ) > 0 [58]. The logarithms L i need to cancel in the final result for the anomalous dimension matrix.

JHEP10(2018)201
β IJ in space-like kinematics. For space-like kinematics of heavy particles (e.g. a top quark incoming and a top quark outgoing), the cusp angle β IJ is real. It can still, however, be chosen to be positive or negative, as the function arccosh has two branches, even for real argument.
The choice of the positive β IJ has at least two advantages: (i) functions that appear in the O α s 2 contribution in eq. (C.11), which we shall discuss in appendix C, are real (in the space-like case) and do not require analytic continuation, as 0 < x < 1 for Li 2,3 (x) and for ln(x), (ii), coth(β IJ ) has a physical interpretation of an inverse of a relative velocity between particles I and J, v IJ , which reads (B.8) The choice β IJ > 0 implies v IJ > 0 and this leads to the following relation between the latter and the β t function introduced in eq. (2.11) v IJ = 2β t 1 + β 2 t . (B.9) We notice that We also note that v IJ ∈ [0, 1], both in the space-like and in the time-like case. Eq. (B.7) has a unique solution for the space-like case with β IJ > 0 which reads and, with help of eq. (B.9) it can be expressed as β IJ in time-like kinematics. We would like to analytically continue β IJ to the time-like region which corresponds to s IJ > 0. The function in eq. (B.13) has a cut at imaginary axis and, depending whether we cross it in the upper or the lower half-plane, we get a different result. The convention is already established in eq. (B.1). It implies that the invariant s IJ has a small positive part which corresponds to analytic continuation over the upper half-plane. And this means that β IJ acquires an imaginary part −iπ in the time-like region. The real part is identical to the one of the space-like kinematics and, at the end, we get T F n f , The QCD β function is given by where the leading coefficient for n f flavours of the active, massless quarks, is The soft renormalization factor Z s can be obtained from the hard renormalization factor Z, determined in ref. [64], by discarding all terms proportional to derivatives of the anomalous dimension and by making the replacement Γ n → γ s n . Up to the order α 2 s , it reads has been introduced in section 2.8, and its explicit result for processes involving massive particles reads [19] γ s qq = C F γ cusp (β 34 , α s ) + 2γ Q (α s ) 1 In order to fully determine the above objects, we need the single-particle massive anomalous dimensions for heavy quarks [58,62]  β coth β + 8C A π 2 6 + ζ 3 + β 2 + coth 2 β Li 3 (e −2β ) + βLi 2 (e −2β ) − ζ 3 + π 2 6 β + β 3 3 + coth β Li 2 (e −2β ) − 2β ln(1 − e −2β ) − π 2 6 (1 + β) − β 2 − β 3 3 , (C.11) and the function g(β) = coth β β 2 + 2β ln(1 − e −2β ) − Li 2 (e −2β ) + π 2 6 − β 2 − π 2 6 . (C.12) The velocity-dependent cusp anomalous dimensions γ cusp (β 34 ) requires careful treatment. It is unambiguously defined for space-like kinematics, in which case the function given in eq. (C.11) is real. In the time-like kinematics, however, γ cusp (β 34 ) has to be analytically continued to the region β ∈ [0, 1], as discussed in appendix B.