1 Introduction

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

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

Guoxing Wang, Xiaofeng Xu, Li Lin Yang and Hua Xing Zhu

School of Physics and State Key Laboratory of Nuclear Physics and Technology,

Peking University, Beijing 100871, China

[0.3cm] Collaborative Innovation Center of Quantum Matter, Beijing, China

[0.3cm] Center for High Energy Physics, Peking University, Beijing 100871, China

[0.3cm] Zhejiang Institute of Modern Physics, Department of Physics, Zhejiang University, Hangzhou 310027, China

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.

## 1 Introduction

Top quark pair production is one of the most important processes at the Large Hadron Collider (LHC). The total cross section at is about . Both the ATLAS and the CMS experiments have collected nearly of integrated luminosity at . This corresponds to 160 million 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 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 next-to-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 production in [30]. Despite the high precisions of these perturbative calculations, the complicated kinematics of 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 . 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 [31, 32]. These include not only the threshold logarithms which arise when the partonic center-of-mass energy approaches the invariant mass , but also the mass logarithms of the form which develop in the boosted region .

The resummation of the threshold logarithms in 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 next-to-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 NLL. 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 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 [46, 47, 34, 35, 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 NLO [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 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.

## 2 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

 Si(x)=Pexp(igs∫0−∞dsvi⋅Aa(x+svi)Tai), (1)

where denotes path ordering, is a 4-vector pointing to the direction of the momentum of the -th parton, which satisfies for massless partons and for massive partons. Note that here we have taken all vectors to be incoming. The boldface is the color generator associated with the -th parton in the color-space formalism [52, 53]. It is evident that the Wilson lines are invariant under the rescaling for any , since this change can be compensated by a change of the integration variable . We could employ this freedom to normalize the direction vectors of massive partons to . This has the physical meaning that is the 4-velocity of the -th parton: where is the mass of the -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 -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

 (2)

where denotes the collection of the directional vectors , is a time-like vector, and 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 [43, 56, 47, 34, 35].

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 [59, 60, 33] for production, in [61, 62] for stop pair production, and in [63, 64] for production. Schematically, we are considering partonic processes of the form

 h1(p1)+h2(p2)→h3(p3)+⋯+hn(pn)+X+Xs(ps), (3)

where and are two incoming massless partons, () are outgoing massive partons, denotes colorless particles such as the Higgs boson and electroweak gauge bosons, and represents additional soft radiations which we want to describe. Here, the momenta () are chosen to be outgoing, but we still keep 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 , which is simply the momentum-space version of the vacuum matrix element of the Wilson loop operator:

 S(ω,{v}) ≡1√d1d2∫dx04πeiωx0/2[W(x,{v})]xμ=(x0,0,0,0) (4)

where the reference vector , 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 and are the dimensions of the representations to which the partons and belong. For later convenience, it is useful to perform a Mellin or Laplace transform into the moment space

 ~s(Λ,{v}) =∫∞0dωexp(−ωΛeγE)S(ω,{v})=1√d1d2W(x0=−2iΛeγE,{v}), (5)

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 . These UV divergences are removed by the operator renormalization

 ~s(L,{v},μ)=Z†s(L,{v},μ)~sbare(Λ,{v})Zs(L,{v},μ), (6)

where is the renormalization scale and . As indicated in the above formula, both the renormalized soft function and the renormalization factor are -dependent and satisfy renormalization group equations (RGEs)

 ddμZs(L,{v},μ) =−Zs(L,{v},μ)Γs(L,{v},μ), (7) ddμ~s(L,{v},μ) =−Γ†s(L,{v},μ)~s(L,{v},μ)−~s(L,{v},μ)Γs(L,{v},μ). (8)

The generic form of the soft anomalous dimension matrix up to the two-loop order can be written as [47]

 Γs(L,{v},μ) =T21+T222γcusp(αs)Ls−∑(I,J)TI⋅TJ2γcusp(βIJ,αs)+∑iγis(αs)+∑IγI(αs) +∑I⎡⎢ ⎢⎣TI⋅T1γcusp(αs)lnv1⋅v2√v2Iv0⋅v2wI1+TI⋅T2γcusp(αs)lnv1⋅v2√v2Iv0⋅v1wI2⎤⎥ ⎥⎦ +∑(I,J,K)ifabcTaITbJTcKF1(βIJ,βJK,βKI) +∑(I,J)∑kifabcTaITbJTckf2⎛⎜ ⎜⎝βIJ,lnwJk√v2IwIk√v2J⎞⎟ ⎟⎠+O(α3s), (9)

where the lower-case indices (, , ) run over massless partons 1 and 2, while the capital indices (, , ) run over massive partons. We have introduced the abbreviations and . The notation denotes unordered tuples of distinct parton indices. The functions and are the famous cusp anomalous dimensions for light-like Wilson lines and time-like Wilson lines, respectively [54, 65, 66]

 γcusp(αs) =αsπ+(αs4π)2[(2689−4π23)CA−809TFNl]+O(α3s), (10) γcusp(b,αs) =γcusp(αs)bcoth(b) +CA2(αsπ)2{π26+ζ3+b2+coth2(b)[Li3(e−2b)+bLi2(e−2b)−ζ3+π26b+b33] (11) +coth(b)[Li2(e−2b)−2bln(1−e−2b)−π26(1+b)−b2−b33]}+O(α3s).

with the cusp angle given by

 βIJ=arccosh(−vI⋅vJ√v2Iv2J−i0). (12)

The single parton soft anomalous dimensions and are

 γis(αs) =(αs4π)2Ci[(−40427+11π218+14ζ3)CA+(11227−2π29)TFNl]+O(α3s), (13) γI(αs) =−CIαs2π+(αs4π)2CI[(−989+2π23−4ζ3)CA+409TFNl]+O(α3s), (14)

where for the fundamental representation, and for the adjoint representation of the gauge group. The three-parton correlation functions and were calculated in [34, 35]. The function describes correlations among three massive partons, and can be written as

 F1(βIJ,βJK,βJI) =(αs4π)243∑L,M,NϵLMNg(βLM)βNLcoth(βNL)+O(α3s), (15)

where the indices (,,) run over (,,) with if (,,) is an even permutation of (,,), and

 g(b)=coth(b)[b2+2bln(1−e−2b)−Li2(e−2b)+π26]−b2−π26. (16)

The function describes correlations among two massive partons and one massless parton, and is given by

 f2⎛⎜ ⎜⎝βIJ,lnwJk√v2IwIk√v2J⎞⎟ ⎟⎠ =−(αs4π)24g(βIJ)×lnwJk√v2IwIk√v2J+O(α3s), (17)

Given the anomalous dimension matrix , one can solve the RGE (7) to obtain the renormalization factor . To this end, it is useful to decompose in (2) into the form

 Γs(L,{v},μ)≡αs4π(A0Ls+Γ0)+(αs4π)2(A1Ls+Γ1), (18)

and then

 lnZs(L,{v},μ) =αs4π(−A02ϵ2+A0Ls+Γ02ϵ) +(αs4π)2[3A0β08ϵ3+−A1−2β0(A0Ls+Γ0)8ϵ2+A1Ls+Γ14ϵ]+O(α3s), (19)

where .

## 3 The soft function for t¯t 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 becomes large. In this paper, we begin with the special case of production, where the partonic processes can be described as

 h1(p1)+h2(p2)→t(p3)+¯t(p4)+Xs(ps), (20)

where with the mass of the top quark. In the soft limit , there are 3 independent Lorentz invariant kinematic variables, which can be chosen as and

 M2≡(p1+p2)2,t1≡(p1−p3)2−m2t. (21)

It is convenient to introduce dimensionless quantities and , defined as

 β=√1−4m2tM2,t1=−M22(1−βy), (22)

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 through the following combinations of the directional vectors :

 v3⋅v1v2⋅v0√v23v1⋅v2 =v4⋅v2v1⋅v0√v24v1⋅v2=−1−βy√1−β2,v3⋅v4√v23v24=1+β21−β2, v3⋅v2v1⋅v0√v23v1⋅v2 =v4⋅v1v2⋅v0√v24v1⋅v2=−1+βy√1−β2, (23)

and we also have

 β34≡arccosh⎛⎜ ⎜⎝−v3⋅v4√v23v24−i0⎞⎟ ⎟⎠=ln1+β1−β−iπ. (24)

To calculate the bare soft function, it is convenient to start from the momentum-space version introduced in Eq. (4). Here we have expressed the dependence on the directional vectors through the quantities and . The perturbative expansion of the momentum-space soft function can be written as

 Sbare(ω,β,y)=δ(ω)1√d1d2+αs4πS(1)bare(ω,β,y)+(αs4π)2S(2)bare(ω,β,y)+⋯, (25)

where denotes the identity operator in color space, and the NLO soft function 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 production in [63, 64], and can be extended to more complicated processes such the simultaneous production of two 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 . It is easy to see that when 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

 (26)

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

 Iij(ϵ,ω,β,y)=∫[dk]vi⋅vjvi⋅kvj⋅kδ(ω−v0⋅k), (27)

where . From symmetry considerations, it is obvious that we only need to calculate , , and . At this point, we note that while the soft function itself does not depend on the normalizations of the directional vectors , it is convenient to fix them in practical calculations. Therefore we will choose the normalizations , and in the following. Note that the normalizations of and are unconventional. In the center-of-mass frame of the incoming partons, these vectors are then parameterized by

 v1=(1,0,0,1),v3=−(1,0,βsinθ,βy), v2=(1,0,0,−1),v4=−(1,0,−βsinθ,−βy). (28)

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

 δ+(k2)≡δ(k2)θ(k0)=12πi(1k2+i0−1k2−i0), (29)

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

 δ(ω−v0⋅k)=12πi(1ω−v0⋅k+i0−1ω−v0⋅k−i0) (30)

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 , and , where is defined as

 Fa1,a2≡∫[dk]δ(ω−v0⋅k)1(v1⋅k)a1(−v3⋅k)a2. (31)

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 with respect to will lead to integrals with 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

 →f(ϵ,β,y,ω)≡(F0,0,F0,1,F1,1)T. (32)

The partial derivative of with respect to then has the form

 ∂β→f(ϵ,β,y,ω)=^A(ϵ,β,y,ω)→f(ϵ,β,y,ω), (33)

where is a matrix. The matrix is a rather complicated function of , , and , which makes the differential equation not so straightforward to solve. It is possible to simplify the above equation by a linear transformation , where the matrix is given by

 ^T(ϵ,β,y,ω)=2Γ(1−2ϵ)π1−ϵω1−2ϵΓ(1−ϵ)⎛⎜⎝1−2ϵ000ϵωβ000ϵω2(1−βy)⎞⎟⎠. (34)

The new vector (so-called “canonical basis”) satisfies a simpler differential equation [76]

 ∂β→g(ϵ,β,y)=ϵ^B(β,y)→g(ϵ,β,y). (35)

Note that the matrix does not depend on and anymore, and is given by

 ^B(β,y)=−^aβ−1+^bβ+^cβ+1−^dβ−1/y, (36)

where

 ^a=⎛⎜⎝000110220⎞⎟⎠,^b=⎛⎜⎝000020402⎞⎟⎠,^c=⎛⎜⎝0001−10−220⎞⎟⎠,^d=⎛⎜⎝000000002⎞⎟⎠. (37)

The vector has no singularity in and can be Taylor-expanded in the form

 →g(ϵ,β,y)=∞∑n=0→g(n)(β,y)ϵn. (38)

The differential equation (35) can then be solved order by order in :

 →g(n+1)(β,y)=∫ββ0dβ′^B(β′,y)→g(n)(β′,y)+→g(n+1)(β0,y), (39)

where the boundary conditions at some boundary point 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

 G(a1,…,an;β)≡∫β0dβ′β′−a1G(a2,…,an;β′), (40)

with . The special case where all the ’s are zero is defined as

 G(0,…,0;β)≡1n!lognβ. (41)

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 at the boundary point , serving as the boundary condition. In our case, it is convenient to choose the point as the boundary. The calculation of the boundary condition can be simplified by observing that the matrix contains a singular term proportional to . It is clear that can produce a coefficient only if itself develops a power-like or logarithmic divergence when . One can easily check that all the master integrals in are regular in the limit . The same applies to the components of since the transformation matrix is also regular. It follows that

 0=limβ→0β∂β→g(ϵ,β,y)=limβ→0βϵ^B(β,y)→g(ϵ,β,y), (42)

which leads to the conditions

 g3(β=0)=−2g1(β=0),g2(β=0)=0, (43)

with being the -th component of . We now only need to directly evaluate the component (which actually doesn’t depend on ) at the boundary , which is very simple:

 g1(ϵ,0,y)=2Γ(2−2ϵ)π1−ϵω1−2ϵΓ(1−ϵ)∫[dk]δ(ω−v0⋅k)=1. (44)

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 . This abstract form is generic and especially useful if we want to apply our result to more complicated processes. However, for practical computations of cross sections, it is convenient to choose a color basis and express the soft function as a matrix in the quark-anti-quark annihilation channel, and a matrix in the gluon fusion channel. Such matrix elements are defined as

 (45)

where is an orthogonal color basis. In accordance with [33], we choose the singlet-octet basis with

 (cq¯q1){a}=δa1a2δa3a4,(cq¯q2){a}=tca2a1tca3a4, (cgg1){a}=δa1a2δa3a4,(cgg2){a}=ifa1a2ctca3a4,(cgg3){a}=da1a2ctca3a4, (46)

where is the color index of the -th parton. We have compared the resulting NLO matrices with those (up to order ) in [33] and find complete agreement.

## 4 The NNLO bare soft function

We now turn to the calculation of the NNLO bare soft function, which is the main new result of our paper. The contributions to the bare NNLO soft function consist of two parts: the virtual-real diagrams and the double-real diagrams. The two-loop virtual diagrams leading to vanishing integrals in dimensional regularization and we do not need to consider them.

### 4.1 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

 ∫[dk1][dk2]δ(ω−v0⋅(k1+k2))F({v},k1,k2), (47)

where denotes the integrand consisting of scalar products among the directional vectors and the two momenta and . 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.

#### One- or two-Wilson-line diagrams

The Feynman diagrams involving only one Wilson line along the vector are depicted in Figure 1. It is clear that such diagrams must be proportional to , and therefore vanish if is light-like. The shaded-blob in the first diagram denotes loops of quarks, gluons and ghosts (we work in the Feynman gauge). Similarly, Figure 2 shows the diagrams involving two Wilson lines in the directions and . The integrals coming from these diagrams can be classified into two integral families. The first family is defined by the following set of 6 propagators:

 {(k1+k2)2,v1⋅k2,v1⋅(k1+k2),v2⋅k1,v3⋅k1,v3⋅(k1+k2)}. (48)

The corresponding integrals have the form

 F(1)a1,a2,a3,a4,a5,a6≡∫[dk1][dk2]δ(ω−v0⋅(k1+k2))6∏i=1(Di)−ai, (49)

where refers to the propagators in Eq. (48). We again feed all the integrals into Reduze2 and FIRE5, which reduce them to 14 master integrals. We collect these master integrals into a vector

 →f(1)(ϵ,β,y,ω)≡( F(1)0,0,0,0,0,0,F(1)1,0,0,0,−1,2,F(1)0,0,0,0,1,0,F(1)0,0,0,0,1,1,F(1)0,1,0,0,0,2, F(1)0,0,1,0,0,2,F(1)0,1,0,0,1,1,F(1)1,1,0,0,1,−1,F(1)1,1,−1,0,1,0,F(1)1,1,0,0,1,0, F(1)1,0,1,0,1,0,F(1)0,0,1,0,1,0,F(1)0,1,1,0,1,1,F(1)1,1,0,1,0,0)T. (50)

We transform the above master integrals to a “canonical basis” via a linear transformation , where the transformation matrix is diagonal with its diagonal entries given by

 ^T(1)(ϵ,β,y,ω)= 8Γ(1−4ϵ)π2−2ϵω3−4ϵΓ2(1−ϵ) ×Diag{ (1−2ϵ)(1−4ϵ)(3−4ϵ),ϵ2(1−2ϵ)ω3β,ϵ(1−2ϵ)(1−4ϵ)ωβ, ϵ2(1−4ϵ)ω2β2,ϵ2ω3β2,ϵ(1−2ϵ)ω3β2,ϵ3ω3β(1−βy),ϵ3ω3, ϵ3ω3β,ϵ3ω4(1−βy),ϵ3ω4(1−βy),ϵ2(1−4ϵ)ω2(1+βy), ϵ3ω4(1−βy)2,ϵ3ω4}. (51)

The transformed vector of master integrals satisfies the differential equation

 ∂β→g(1)(ϵ,β,y)=ϵ⎛⎝−^a(1)β−1+^b(1)β+^c(1)β+1+^d(1)β+1/y−^e(1)β−1/y⎞⎠→g(1)(ϵ,β,y), (52)

where ,