Master integrals for the NNLO virtual corrections to $q \bar{q} \rightarrow t \bar{t}$ scattering in QCD: the non-planar graphs

We complete the analytic evaluation of the master integrals for the two-loop non-planar box diagrams contributing to the top-pair production in the quark-initiated channel, at next-to-next-to-leading order in QCD. The integrals are determined from their differential equations, which are cast into a canonical form using the Magnus exponential. The analytic expressions of the Laurent series coefficients of the integrals are expressed as combinations of generalized polylogarithms, which we validate with several numerical checks. We discuss the analytic continuation of the planar and the non-planar master integrals, which contribute to $q {\bar q} \to t {\bar t}$ in QCD, as well as to the companion QED scattering processes $ e e \to \mu \mu$ and $e \mu \to e \mu$.


Introduction
The top quark is the heaviest elementary matter particle. Its large production rate at the CERN LHC enables precision studies, thereby testing the Standard Model of particle physics at an unprecedented level, and potentially uncovering indirect evidence for new physics effects. Precision measurements of top quark observables [1-4] must be confronted with equally precise theory predictions, thereby requiring the perturbation theory descriptions of these observables to be extended to high orders.
The available second-order (next-to-next-to-leading order, NNLO) QCD corrections to the top quark pair production process (initially for the total cross section [5] and subsequently for differential distributions [6][7][8][9]) deliver a competitive level of theory predictions, and are enabling a multitude of precision studies with top quark pairs. A key ingredient to these calculations are the two-loop QCD corrections to the matrix elements for top pair production in quark-antiquark annihilation and gluon fusion. While numerical representations of these two-loop matrix elements were derived already some time ago [10,11], only partial analytical results are available for them up to now [12][13][14][15] (and used in partial validations [16,17] of the total NNLO cross section calculation). Such analytical results allow in-depth investigations into the structure of the matrix elements, enabling to investigate limiting behaviors and analyticity structures, as well as providing an independent approach to their numerical evaluation. Up to now, full results for the two-loop top quark -1 -production matrix elements could not be derived due to incomplete knowledge on the relevant two-loop Feynman integrals. Important steps towards the evaluation of the missing two-loop functions have been recently undertaken in [18][19][20][21].
In this work, we complete the task of determining all the non-planar two-loop functions that are needed for the evaluation of the scattering amplitudes for the quark initiated channel qq → tt at NNLO in QCD. Owing to the large value of the top mass compared to the mass of the two incoming quarks, which we suppose to have different flavor, we treat the latter as massless.
The results of this paper represent an additional milestone of the research program dedicated to the two-loop QCD and QED corrections to the interaction of two fermionic currents, that was initiated with the study of the muon-electron scattering, in the context of the MUonE experiment [22,23], which is currently under evaluation at CERN. By following the same the path of the calculation of two-loop integrals for µe → µe and the crossing-related processes considered in [24,25], we adopt a consolidated strategy [26,27], which was proven to be particularly effective in the context of multi-loop integrals that involve multiple kinematic scales [24,25,[27][28][29][30][31]. By means of integration-by-parts identities (IBPs) [32][33][34], we identify a set of 52 master integrals (MIs) that we evaluate analytically, through the differential equations method [35][36][37][38]. In particular, we consider an initial set of MIs that obey a system of first-order differential equations (DEQs) in two independent kinematic variables that is linear in the space-time dimensions d. The system is subsequently cast in canonical form [26] by means of the Magnus exponential matrix [27]. The matrix associated to the canonical system, where the dependence on (d−4) is factorized from the kinematics, is a logarithmic differential form, which -once parametrized in terms of suitable variables -has a polynomial alphabet, constituted of eleven letters. Therefore, the canonical MIs are found to admit a Taylor series representation around d = 4, whose coefficients are combinations of generalized polylogarithms (GPLs) [39][40][41][42]. The otherwise unknown integration constants are determined either from the knowledge of the analytic expression of the MIs in special kinematic configurations or by imposing their regularity at the pseudo-thresholds of the DEQs. Finally, we show how the MIs, that we initially compute in an unphysical region, can be analytically continued to the top-pair production region. Due to the non-trivial structure of their branch-cuts, the analytic continuation of the two-loop functions considered in this paper represents a paradigmatic case, that can be useful for the study of other planar and non-planar diagrams that involve massive particles. As a byproduct of the current analysis, we obtain the analytic continuation to the physical region of the functions required for the µe → µe and ee → µµ scattering in QED [24,25] In the completion of this calculation, we benefited from publicly available software dedicated to multi-loop calculus. The IBPs decomposition and the generation of the dimensionshifting identities and DEQs obeyed by the MIs have been performed with the packages Kira [44], LiteRed [45,46] and Reduze [47,48]. The analytic expressions of the MIs, which we evaluate numerically with the help of GiNaC [49], were successfully tested against the  Figure 1: Representative non-planar diagrams contributing at two-loops to q(p 1 )+q(p 2 ) → t(p 3 ) +t(p 4 ) (top row), and the associated integral families (bottom row). Massive propagators and external legs are depicted using thick lines. The diagrams have been drawn with axodraw2 [52].
numerical values provided either by SecDec [50] or, for the most complicated 7-denominator topologies, by an in-house algorithm. Beside these important validations, our results have been successfully compared against the computation of the master integrals relevant to the same integral topologies expressed in a different basis set, independently obtained by Becchetti, Bonciani, Casconi, Ferroglia, Lavacca and von Manteuffel [51], published in tandem to the current manuscript.
The remainder of the paper is organized as follows: in section 2 we set our notation and conventions for the non-planar two loop functions. In section 3, we present the system of DEQs obeyed by the MIs and their solutions in the unphysical region. In section 4, we study the analytic continuation of the MIs to the physical region. Finally, in section 5, we discuss the numerical evaluation of the 7-denominator integrals. Appendix A contains the coefficients of the linear combinations of MIs that satisfy a canonical system of DEQs.
The analytic expressions of the considered MIs are given in the ancillary files accompanying the arXiv version of this publication.

The non-planar four-point topology
In this paper, we complete the determination of the Feynman integrals for the qq → tt scattering process with kinematics specified by where m is the top quark mass. The full set of two-loop three-point integrals that are involved in the process has been known for some time [53][54][55][56][57][58][59][60], as well all the relevant planar four-point functions [12,13,24]. On top of these contributions, the evaluation of the full double-virtual matrix element requires the computation of two-loop non-planar Feynman diagrams. Representative non-planar diagrams are shown in the top row of figure 1. In the bottom row, we also show the integral families onto which we map those diagrams. Massive propagators and external legs are depicted using thick lines. The MIs for the QED-like family A 1 are already available in the literature, as they have been studied in the context of the NNLO QED corrections to eeµµ processes (with suitable redefinitions of the momenta and of the Mandelstam invariants). In particular, they have been first computed in [25] (in an unphysical region, to be analytically continued), and later in [43] (directly in the heavy-fermion-production kinematic region). As for the genuine QCD contributions, the MIs for family N 1 have been computed in [61], while the MIs for family N 2 are the subject of the present publication.
In arbitrary d dimensions, we parametrize the integrals of family N 2 as where D i are the inverse scalar propagators with k 1 and k 2 denoting the loop momenta. The analytic calculation described in section 3 is performed by expanding the MIs around d = 4, while the numerical evaluation presented in section 5 as a check of our results is carried over around d = 6. Therefore, we set ≡ (d * − d)/2, where d * = 4 and d * = 6 according to the case. In addition, we define our integration measure as where µ is the 't Hooft scale of dimensional regularization. In this convention, the two-loop tadpole integral 2 I [4−2 ] (0, 2, 0, 0, 0, 2, 0, 0, 0) is normalized to 1.
3 Solution of the system of differential equations where the constraint on the Mandelstam invariants s + t + u = 2m 2 is understood. The dimensionless variables w and z appearing in eq. (3.1) were already introduced in [25] for the computation of the MIs of the non-planar family A 1 (modulo the relabeling s ↔ t). Also in this case, the above change of variables rationalizes the coefficients of the canonical DEQs, hence allowing to express the MIs in terms of GPLs.

Differential equations for master integrals
We determine a canonical basis of MIs in d = 4 − 2 through the Magnus exponential algorithm described in [27,28]. As a starting point, we identify a set of 52 independent integrals F i , whose DEQs depend linearly on the dimensional regularization parameter , where the T i are the integrals graphically represented in figures 2 and 3.
In particular, the numerators of integrals F 49,...,52 , are found by following the ideas in [62], i.e. by determining a set of canonical integrals through the inspection of their Figure  four-dimensional maximal-cuts. To this aim, we first localize the integral over k 2 , which corresponds to the non-planar part of the diagram specified by D 2,3,6,7 . By enforcing the additional constraints k 2 1 =0 and s = 2k 1 · (p 3 + p 4 ) (which originate from the cut of the Figure  propagators depending on k 1 ), we obtain where we have denoted δ i = δ(D i ) and assumed the integral to contain some arbitrary numerator depending on k 1 . From eq. (3.3), we observe that the maximal-cut of the oneloop subdiagram exposes two hidden propagators, The latter, together with the residual uncut propagators D 1,4,5 , form a one-loop pentagon integral, known to obey non-canonical DEQs around four dimensions. Therefore, we choose the numerator factor N (k 1 ) such that it cancels one or both the hidden propagators, resulting in either a box or triangle integral, which both satisfy canonical DEQ. In this way, we determine three out of the four numerators corresponding to the integrals F 49,...,51 , as they are displayed in figure 3. The last numerator F 52 is defined to contain, besides the factor D 10 , also the auxiliary denominator D 9 , which, depending on k 2 , ensures the linear independence from the other three basis integral of the sector, without spoiling the properties of the DEQs.
Once a basis of MIs with -linear DEQs has been determined, the integrals of eq. (3.2) can be rotated into a canonical basis of MIs I i by means of the Magnus exponential matrix. Through this procedure, we find that the integrals obey canonical DEQs in both kinematic variables w and z. In eq. (3.4) we have used the compact notation λ s = √ −s √ 4m 2 − s. The analytic expression of the coefficients c i,52 , that are here omitted for brevity, can be found in appendix A, as well as in the ancillary files of the arXiv version of the manuscript, which contain eq. (3.2) and eq. (3.4) in electronic format.
By organizing the 52 MIs into a vector I( , w, z) and by combining the DEQs in w and z obeyed by the latter into a single total differential, we can write the canonical DEQs compactly as dI = dAI , where dA is the 1-form with M i being constant matrices (that are provided as ancillary files in the arXiv submission of this paper). The alphabet of the DEQs, i.e. the set of arguments η i of the d log-form in eq. (3.6), can be taken to be formed by the same 12 letters that appear in the calculation of the MIs for the QED-like topology A 1 presented in ref. [25]: For the MIs of the topology N 2 under consideration, the matrix M 11 vanishes identically. Nevertheless, we adopted the above notation for consistency with ref. [25]. The analytic expression of the MIs is first derived in the region of the wz-plane where the whole alphabet is real and positive, -9 -which corresponds to the unphysical kinematic region The analytic continuation to the full kinematic plane is discussed in section 4. By definition, the integrals given in eq. (3.4) are finite in the → 0 limit. Therefore, the vector I( , w, z) can be expanded as a Taylor series around = 0 as (3.10) The n-th order coefficient of eq. (3.10) is given by where the I (i) (w 0 , z 0 ) are constant vectors and ∆ (k) is the operator that represents k path-ordered iterated integrations of dA along a piecewise-smooth path γ in the wz-plane. As the roots of the rational alphabet given in eq. (3.7) are purely algebraic, the iterated integrals of eq. (3.12) can be directly mapped into combinations of GPLs, G( a n ; x) ≡ G(a 1 , a n−1 ; x) ≡ x 0 dt 1 t − a 1 G( a n−1 ; t), (3.13) (3.14) The length n of the vector a n is commonly referred to as the transcendental weight of G( a n ; x) and it amounts to the number of repeated integrations defining the GPL. We found it convenient to determine the solution of the DEQs in terms of GPLs, up to O( 5 ) terms, by integrating first in w and then in z. Consequently, the GPLs that appear in the analytic expression of the MIs fall into two classes: GPLs with argument w and weights drawn from the set and GPLs with argument z and with weights drawn from the set {0 , ±1}. Due to the positivity of the alphabet, the combinations of GPLs that appear in the solution I( , w, z) are real-valued in the region defined by eq. (3.8). It therefore follows that all the imaginary parts associated to the MIs whose physical thresholds are overstepped in this region are made explicit in the integration constants I (i) (w 0 , z 0 ).

Boundary conditions
In order to completely specify the analytic expression of the MIs, a suitable set of boundary conditions must be imposed on the general solution of the system of DEQs. Boundary conditions are determined by imposing the regularity of the MIs at the pseudo-thresholds of the DEQs, that entails, order by order in expansion, a linear relation between the MIs. Such regularity conditions are complemented by the knowledge of the analytic expression of a limited number of input integrals in special kinematic configurations. In the case under study, the boundary constants of all the MIs can be expressed as combinations of GPLs of argument 1 and complex weights, the latter arising from the specific kinematic limits imposed on the alphabet of eq. (3.7). With the help of GiNaC, we were able to numerically verify that, at each order in , such combinations of constant GPLs are proportional to uniform combinations of the transcendental constants π, ζ k and log 2.
The boundary constants of each integral have been determined as follows: • the integrals I 1,...,7,10,...,14,32,...,37 are either common to the two-loop topologies discussed in reference [24,25] (to which we refer the reader for the discussion of the boundary value fixing procedure) or related to them by simple kinematic crossing, i.e. by some interchange of the Mandelstam invariants; • the boundary constants of I 8, 9,18 have been fixed by demanding finiteness of the MIs in the limit s → 0; • the boundary constants of I • the boundary constant of I 45 has been fixed by demanding the finiteness of the integral • the boundary constant of I 48 has been fixed by taking the massless limit, as described in [24].
We provide the analytic expressions of the MIs in electronic form in the ancillary files attached to the arXiv submission of the paper.
-11 - The results of section 3 have been obtained in the unphysical region s, t < 0. Therefore, the analytic continuation of such expressions to the tt production kinematics must be performed. In terms of the Mandelstam invariants, this region is defined by where the boundaries of the allowed interval for t are in one-to-one correspondence, in the center-of-mass frame, with the minimum and maximum scattering angles of the tt pair with respect to the beam line. For completeness, we also quote the physical region for elastic scattering, corresponding to the crossed t-channel process: In the case of non-planar four-point integrals, the analytic continuation of the GPLs is quite non trivial. As originally noted in [63,64], thresholds associated with all the Mandelstam invariants appear simultaneously, and s, t, u should be treated as independent variables when discussing the analytic continuation of the expressions for the MIs. On the other hand, the approach we follow enforces the constraint s + t + u = 2m 2 from the outset, yielding a system of DEQs in two variables, e.g. s and t. One way out could be to enforce the Mandelstam constraint only at a later stage (see e.g. [65]), at the price of considerably complicating the problem to be solved due to the presence of an extra scale.
In this paper we address the analytic continuation in a different way by exploiting the iterated-path-integral nature of our canonical MIs, together with the so-called firstentry condition [66,67], in order to devise an effective prescription. Our approach allows to analytically continue the MIs everywhere in the kinematic plane, and in particular to evaluate our results in the tt production region. As a byproduct of our analysis, we also obtain the analytic continuation of the MIs for µe scattering, presented in [25], both to the di-muon production region and to the elastic scattering region.
We already observed in section 3 that our canonical basis of MIs can be expressed, order by order in , as a linear combination of GPLs and constants of uniform weight. From the first-entry condition it follows that only the innermost integration contributes to the discontinuities of the MIs. Strong restrictions on the analytic structure of the MIs are imposed already at the level of the canonical DEQs (by the coefficient matrices in the d log form), but knowledge of the boundary conditions is essential to fully pin it down. By inspection of our result (computed up to weight 4 in the region s, t < 0), we find that the GPLs originating from the innermost integration are the following: G 0 (z), G 0 (w), G 1 (w), G z 2 (w). This can be traced back to the fact that, of all the letters η i that appear in dA (see eq. (3.7)), only four contribute to the first integration, namely η 1, 3,4,9 . Quite remarkably, one can build four combinations of the η 1,3,4,9 that correspond to simple functions of the Mandelstam invariants, whose logarithms exhibit the branch cuts expected from the normal thresholds -12 -of the four sunrise sub-topologies. If we define where use has been made of the relation s + t + u = 2m 2 , then one can easily find the following GPL representations for the corresponding logarithmic differentials We refer to log θ 1,2,3,4 as physical logarithms. In figure 4 we show the physical regions for the s-channel production and the t-channel scattering processes, together with the region in which we solved the system of differential equations, and the thresholds of the physical logarithms. For completeness we also give a more transparent rearrangement of the other letters: η 10 /η 9 ≡ θ 6 = −t/m 2 , (4.12) η 2 2 /η 1 ≡ θ 7 = 4 − s/m 2 , (4.13) η 2 3 η 5 η 6 η 7 η 8 /(η 1 η 2 9 ) ≡ θ 8 = 1 − tu/m 4 , (4.14) One can prove that, in the region s, t < 0, eqs. (4.7)-(4.10) also hold if the total differential operator is dropped, without adding any integration constants. By choosing a suitable analytic continuation prescription on the Mandelstam invariants one can then evaluate the integrated expressions in the full kinematic plane, in an unambiguous way. One can then check whether those expressions reproduce the imaginary parts of the corresponding physical logarithms. The simple prescription we adopt is defined in two steps:  Figure 4: In the plot we show a representative portion of the kinematic phase space, parametrized in terms of (s, t). The region in which we solved the system of differential equations, s, t < 0, is marked in green. The physical region for the s-channel production process, given in eq. (4.1), is highlighted in blue. In orange we also show the physical region for the t-channel process, given in eq. (4.2), which is relevant for µe scattering. The dashed lines indicate the thresholds of the physical logarithms in eq. (4.7)-(4.10).
1. As for the Mandelstam invariants, we express the real part of u in terms of the real parts of s and t, for which we use the standard Feynman prescription, but we give u an independent prescription (i.e. before using u(s, t) = 2m 2 − s − t in the definition of z, eq. (3.1)) where iε is an infinitesimal positive imaginary part, θ(x) is the Heaviside step function, and the presence of the constant iε term in the last equation guarantees that the evaluation of the GPLs on spurious branch cuts (that are developed even for s, t < 0) -14 -is always unambiguous. It can be shown, by repeated application of the identity log ab = log a + log b + 2πi [θ(− Im a)θ(− Im b)θ(Im ab) −θ(Im a)θ(Im b)θ(− Im ab)] , (4.21) that the above prescription is sufficient to reproduce, with the GPL representation in the variables (w, z) of eqs (4.7)-(4.10), the physical logarithms log θ 1,2,3 (that, for instance, completely determine the analytic structure of the s-channel sunrise MIs and the t-channel sunrise MIs).
2. It remains to be verified whether the above prescription allows to correctly reproduce also the imaginary part of log θ 4 (the one that, for instance, determines the discontinuity of the u-channel sunrise MIs across the one-particle branch cut). This is not guaranteed since, as stressed at the beginning of this section, we only have two independent variables at our disposal, while having to deal simultaneously with thresholds in all the three channels. The virtue of our prescription (4.18)-(4.20) is that the representation of log θ 4 in terms of GPLs (corresponding to eq. (4.10)) • for u > m 2 is always on the physical Riemann sheet; • for u < m 2 always ends up on the wrong side of the branch cut, i.e. the imaginary part is always −iπ instead of +iπ.
We can therefore apply, as a second step, a simple correction to the combination of GPLs corresponding to log θ 4 (and not to the other three combinations) to bring our GPL representation for the MIs to the correct Riemann sheet. At weight one one can perform the replacement which is then propagated iteratively to higher weights.
All the explicit imaginary constants in our solution, as stated in section 3, originate from the (iterated) integrations over d log θ 4 . Indeed we integrate our DEQs in the region where s, t < 0, so that u > 2m 2 . The combinations of GPLs corresponding to such integrations (at any weight) are then always accompanied by a constant term, namely an additional −iπ. Therefore, the net effect of the correction (4.22) is to flip the sign of the imaginary constants in our solution. In summary, our effective way of implementing the analytic continuation of the result for the full set of MIs is of explicitly [64], the first-entry condition guarantees that the analytic continuation of the full set of MIs at all weights is also correctly obtained. In particular, it is not necessary to make sure the GPL representation also reproduces the imaginary parts coming from the evaluation close to the branch cuts of the logarithms of the unphysical letters, eqs. (4.11)-(4.17). It is instead sufficient to always introduce an "auxiliary" iε prescription in order to avoid the ambiguous evaluation on such branch cuts (in our case this is inherited from (4.18)-(4.20)). Our strategy for the analytic continuation has been validated by thorough numerical checks performed either with the help of SecDec or with the techniques outlined in section 5.
It is now clear that we can also obtain, by the same argument, the analytic continuation of the results for the MIs of the QED-like topology A 1 (see figure 1) presented in [25]. The only difference with respect to the present case (besides a trivial relabeling of the Mandelstam invariants s ↔ t to match the notations), is that the letter η 11 contributes to the d log form with a non-zero coefficient matrix (but never appears in the first entry), and that η 1 = θ 1 is not a physical logarithm anymore (as expected due to the absence of a two-massive-particles normal threshold). Since the analytic continuations of log θ 1 and log θ 3 are not independent, in practice this difference does not change the situation, and the same procedure described above can then be used for the analytic continuation of the MIs from s, t < 0 to the physical regions for the µe → µe and e + e − → µ + µ − processes, as confirmed also in this case by our numerical checks.
We stress that the method outlined above is fully general, as it only relies on the analyticity properties of the canonical MIs, and on their iterated-integral representation. It is in particular independent of the presence of massive propagators or massive external legs.

Numerical validation of the non-planar box integrals
Using the analytic continuation as described in the previous section, the expressions for our MIs have been numerically evaluated in several points across the whole phase space, including the unphysical region s, t < 0 and the region relevant for tt production, eq. (4.1). In order to cross-check our analytic calculation, we numerically computed the MIs (or linear combinations of the latter) in some benchmark points with an alternative method, namely by integrating directly their Feynman-parameter representation. In particular, the integrals I i with i = 1, . . . , 48 were computed with the package SecDec. For the most complex topologies, corresponding to the non-planar box integrals I i with i = 49, . . . , 52, we used Reduze to identify an alternative set of independent MIs that are quasi finite [48] in d = 6. On the one hand, the latter have been computed semi-numerically by means of an in-house algorithm [25]. On the other hand, these integrals can be analytically related to our set of MIs by dimension-shifting identities [68,69] and IBPs, implemented in LiteRed, therefore allowing for a numerical comparison.
The definition of the 4 non-planar 7-denominator MIs that are finite in d = 6 dimensions, together with our results at the phase-space point s = −1/7, t = −1/3, m 2 = 1, are collected in table 1. In the next subsection, we use the first of those integrals as an where we used the notation x i 1 i 2 ...in = x i 1 +x i 2 +. . .+x in and the abbreviation Γ ≡ Γ(1+ ). Now we perform as many analytic integrations as possible. In particular, we eliminate the δ-function by integrating over x 3 , and we make the change of variables In this way, the polynomial ∆ becomes linear in x 4 , and eq. (5.4) is rewritten as where The integral over x 4 in eq. (5.6) is finite for d → 6. In this limit, we get 0 dx 1 f 2 − g 2 x 1 ln g 2 x 2 1 + g 1 x 1 + g 0 f 3 x 1 + P 6 , , Finally, we integrate over x 1 and reduce eq. (5.8) to where Differently from the case of the non-planar integral of the family A 1 , which was evaluated with a similar strategy in [25], the integrand contains 6 distinct dilogarithms, whose arguments depend of the algebraic functions P 1 , P 2 , P 3 , P 4 , R 1 and R 2 , which contain the square root of the same polynomial.

Numerical integrations
We rescale the four remaining integration variables in eq. (5.10) in order to map them onto a four-dimensional hypercube of unit side, In this way, the new variables t i have to be integrated over [0, 1]. The six dilogarithms that appear in eq. (5.10) have branch points, which correspond to the hypersurfaces defined by the equations In order to obtain a fast convergence of the numerical integration, the integrands must be sampled carefully in the neighborhood of these branch points. We start from the integration over t 4 . We split the integration interval at the N 4 (t 1 , t 2 , t 3 ) solutions z 4j (t 1 , t 2 , t 3 ) of eq. (5.13), which satisfy 0 ≤ z 4j ≤ 1.
dt 4 , z 40 = 0 , z 4N 4 = 1 . (5.14) Now we consider the integration over t 3 . We split the integration interval at the N 3 (t 1 , t 2 ) zeros z 3j (t 1 , t 2 ) of the discriminants (polynomials in (t 1 , t 2 , t 3 )) that appear in the zeros z 4j and which satisfy 0 ≤ z 3j ≤ 1. These are the points where the hypersurfaces of eq. (5.13) are tangent to the hyperplane t 4 = constant, Next, we consider the integration over t 2 . We split the integration interval at the N 2 (t 1 ) zeros z 2j (t 1 ) of the discriminants (polynomials in (t 1 , t 2 )) that appear in the zeros z 3j and which satisfy 0 ≤ z 2j ≤ 1, Finally, the integration interval in u is subdivided in 2 n subintervals, and Gauss-Legendre integration over 16 points in each subinterval is employed. All the singularities in the integrands are integrable logarithmic singularities. Therefore, we can safely set a very small value of ε, like 10 −30 , adequate for 16-digit calculations. By using 16 subdivisions in each interval and in every variable, we find that our integral, in the phase space point s = −1/7, t = −1/3, m 2 = 1, has the value We adopt a similar procedure for the other integrals shown in table 1. In all cases, after the analytic integrations the integrands, in the d → 6 limit, are found to contain combinations of logarithms of ratios of the polynomials P i . Henceforth, the decomposition of the integration domain, as well as the numerical integration are carried out along the same lines described above as for the scalar box integral.
-20 -In this paper, we presented the analytic expressions of the master integrals for a set of non-planar two-loop Feynman graphs, with two quark currents exchanging massless gauge bosons. Our results are the last missing ingredients required for the analytic evaluation of the double-virtual contribution to the scattering amplitude for the process qq → tt at NNLO in QCD, which was so far known only numerically. The present computation completes the calculation of all the required master integrals, hence proving that the analytic evaluation of such amplitude is indeed feasible.
The two-loop integrals were evaluated by means of the differential equations method, which, combined with the ideas of the Magnus exponential matrix and of the canonical basis, yielded a representation of the master integrals in terms of generalized polylogarithms. The canonical systems of differential equations was conveniently solved in a non-physical region. Subsequently, we studied, in the presence of massive internal lines and of a nontrivial structure of branch cuts, the analytic continuation of the two-loop functions to the physical region relevant for the process under consideration.
The results of this paper represent an important milestone of the research program dedicated to the evaluation of integrals originating from the planar and non-planar twoloop four-point diagrams that contribute to both QED and QCD processes, which include, beside the top-pair production at hadron colliders, also di-muon production at lepton colliders, as well as muon-electron elastic scattering, which is the investigation target of the novel experiment MUonE, recently proposed at CERN.