Master Integrals for the two-loop, non-planar QCD corrections to top-quark pair production in the quark-annihilation channel

TL;DR

This paper analytically computes master integrals for two-loop non-planar QCD corrections in top-quark pair production, advancing precision calculations in particle physics.

Contribution

It provides the first analytic solutions for these complex integrals using differential equations and harmonic polylogarithms, enhancing theoretical predictions.

Findings

Explicit epsilon expansions of integrals up to weight four

Use of generalized harmonic polylogarithms for solutions

Facilitation of precise top-quark pair production calculations

Abstract

We present the analytic calculation of the Master Integrals for the two-loop, non-planar topologies that enter the calculation of the amplitude for top-quark pair hadroproduction in the quark-annihilation channel. Using the method of differential equations, we expand the integrals in powers of the dimensional regulator $\epsilon$ and determine the expansion coefficients in terms of generalized harmonic polylogarithms of two dimensionless variables through to weight four.

Equations

$$s=(p_1+p_2{)}^2,t=(p_1-p_3{)}^2,u=(p_1-p_4{)}^2,$$

$$s>4m^2,t=m^2-\frac{1}{2}\left(s-\sqrt{s(s-4m^2)}\cos \theta \right),$$

$$\int {\mathcal{D}}^d{k}_1{\mathcal{D}}^d{k}_2\frac{D_4^{-a_4}{D}_6^{-a_6}}{D_1^{a_1}{D}_2^{a_2}{D}_3^{a_3}{D}_5^{a_5}{D}_7^{a_7}{D}_8^{a_8}{D}_9^{a_9}},$$

$$\int {\mathcal{D}}^d{k}_1{\mathcal{D}}^d{k}_2\frac{D_5^{-b_5}{D}_6^{-b_6}}{D_1^{b_1}{D}_2^{b_2}{D}_3^{b_3}{D}_4^{b_4}{D}_7^{b_7}{D}_8^{b_8}{D}_9^{b_9}}.$$

$${\mathcal{D}}^d{k}_i=\frac{d^d{k}_i}{i{\pi }^{\frac{d}{2}}}{e}^{ϵ{\gamma }_E}{\left(\frac{m^2}{{\mu }^2}\right)}^{ϵ},$$

$$D_i$$

$${\partial }_i\vec{f}(\vec{x},ϵ)=A_i(\vec{x},ϵ)\vec{f}(\vec{x},ϵ),$$

$${\partial }_i{A}_j-{\partial }_j{A}_i+\left[A_i,A_j\right]=0,$$

$$d\vec{f}(\vec{x},ϵ)=ϵd\tilde{A}(\vec{x})\vec{f}(\vec{x},ϵ).$$

$$\vec{f}(\vec{x},ϵ)=\mathbb{P}\exp \left(ϵ{\int }_{\gamma }d\tilde{A}(\vec{x})\right){\vec{f}}_0(ϵ),$$

$$G({\alpha }_1,\dots ,{\alpha }_n;z)={\int }_0^z\frac{dt}{t-{\alpha }_1}G({\alpha }_2,\dots ,{\alpha }_n;t),$$

$$G({\alpha }_1;z)={\int }_0^z\frac{dt}{t-{\alpha }_1}\text{mbox}for{\alpha }_1\ne 0,\text{mbox}andG({\vec{0}}_n;z)=\frac{{\log }^n(z)}{n!},$$

$${\mathcal{T}}_1^A,{\mathcal{T}}_1^B$$

$$t$$

$${\mathcal{T}}_{36}^A,{\mathcal{T}}_{33}^B$$

$${\mathcal{T}}_{29}^A$$

$$s$$

$${\mathcal{T}}_{26}^A,{\mathcal{T}}_{23}^B$$

$${\mathcal{T}}_{44}^A,{\mathcal{T}}_{38}^B$$

$${\mathcal{T}}_{49}^A$$

$$f_1^A$$

$$f_2^A$$

$$f_3^A$$

$$f_4^A$$

$$f_5^A$$

$$f_6^A$$

$$f_7^A$$

$$f_8^A$$

$$f_9^A$$

$$f_{10}^A$$

$$f_{11}^A$$

$$f_{12}^A$$

$$f_{13}^A$$

$$f_{14}^A$$

$$f_{15}^A$$

$$f_{16}^A$$

$$f_{17}^A$$

Full text

††preprint: MSUHEP-18-021

1 Introduction

Theoretical predictions for top-antitop pair production at hadron colliders are known in perturbative QCD up to next-to-next-to-leading order (NNLO) Baernreuther:2012ws ; Czakon:2012zr ; Czakon:2012pz ; Czakon:2013goa ; Czakon:2014xsa ; Czakon:2015owf ; Czakon:2016dgf ; Czakon:2016ckf . Recently, also the NLO electroweak corrections to this process were evaluated Czakon:2017wor . Predictions at NNLO in QCD are available for the total cross section and for distributions that are differential with respect to quantities which depend on the momenta of the top-antitop pair, such as the pair invariant mass, the top (or antitop) transverse momentum and rapidity, etc.

From the technical point of view, the numerical calculations carried out in Baernreuther:2012ws ; Czakon:2012zr ; Czakon:2012pz ; Czakon:2013goa ; Czakon:2014xsa ; Czakon:2015owf ; Czakon:2016dgf ; Czakon:2016ckf represent a landmark in the field of the evaluation of higher-order corrections in perturbative QCD. One of the main technical problems that was necessary to solve in order to achieve NNLO accuracy was the evaluation of two-loop $2\to 2$ amplitudes with massive and massless propagators. The evaluation had to be carried out for arbitrary values of the Mandelstam invariants $s$ and $t$ and of the top-quark mass $m$. The problem was solved by evaluating numerically these diagrams in a grid of points covering all of the physics phase space in the $s-t$ plane, for a fixed value of $m$. The evaluation in each single point was carried out by solving numerically differential equations satisfied by the Master Integrals (MIs) present in the problem. The numerical solution of large sets of differential equations is not only technically challenging but it also requires a significant amount of CPU time. In addition, it was necessary to evaluate analytically the boundary conditions to be used in the numerical solution of the differential equations.

In this context, an analytic calculation of the two-loop amplitudes contributing to top-quark pair production has a twofold purpose: on the one hand, it provides an independent check of the results obtained numerically; on the other hand it could provide a faster and cheaper (in terms of CPU time) way to evaluate the two-loop corrections needed in order to obtain phenomenological predictions for this process.

A complete analytic computation of the top-pair production cross section to NNLO in QCD is not yet available, although many of the necessary elements were evaluated in the recent past. In particular, the matrix elements for the one-loop $2\to 3$ process are known Dittmaier:2007wz ; Bevilacqua:2010ve ; Bevilacqua:2011aa ; Melnikov:2010iu . Furthermore, progress was also made in the determination of infra-red (IR) subtraction terms which are needed to regularize IR divergences in collinear and soft regions of the phase space during the integration Abelof:2014fza ; Abelof:2014jna ; Abelof:2015lna ; Bonciani:2015sha . Very recently the computation of the NNLO IR subtraction terms was completed in the $Q_{T}$ subtraction formalism Angeles-Martinez:2018mqh ; Catani:2019iny . Finally, the one-loop squared matrix elements were calculated in Korner:2008bn ; Kniehl:2008fd ; Anastasiou:2008vd . Analytic results for the interference between two-loop $2\to 2$ diagrams and tree-level amplitudes are available only in part.

Two-loop contributions to the $t\bar{t}$ production process in hadronic collisions are required for two partonic channels: $q\bar{q}\to t\bar{t}$ (quark-annihilation channel) and $gg\to t\bar{t}$ (gluon fusion channel). The interference of the two-loop amplitude in the quark-annihilation channel with the corresponding tree-level amplitude can be expressed in terms of ten gauge independent functions. Each one of these functions is proportional to a different color coefficient. In the rest of this work we refer to these functions as “color factors”. The color structure in the gluon-fusion channel is more complicated, and it can be expressed in terms of sixteen color factors.

All of the ten color factors in the $q\bar{q}$ channel are known numerically Czakon:2008zk and their infrared poles are known analytically Ferroglia:2009ep ; Ferroglia:2009ii . For eight out of the ten color factors a complete analytic expression, written in terms of generalized harmonic polylogarithms (GPLs) Goncharov:polylog ; Goncharov2007 ; Remiddi:1999ew ; Vollinga:2004sn , was found in Bonciani:2008az ; Bonciani:2009nb . The remaining two color factors in the quark-annihilation channel are not known analytically to date.

All of the sixteen color factors appearing in the two-loop corrections in the gluon-fusion channel are known numerically Baernreuther:2013caa and the analytic expression of all the infrared poles was evaluated in Ferroglia:2009ep ; Ferroglia:2009ii . In addition, a complete analytic expression (again written in terms of GPLs) is known for ten out of the sixteen color factors in the gluon fusion channel Bonciani:2010mn ; vonManteuffel:2013uoa ; Bonciani:2013ywa . The remaining six color factors in this partonic channel are known to involve elliptic integrals. Very recently, the MIs for planar topologies that involve a closed heavy fermionic loop were studied in Adams:2018bsn ; Adams:2018kez . These MIs contribute to one of the six gluon-channel color factors that are not known analytically.

In this paper we focus on the analytic calculation of the MIs that are needed to complete the evaluation of the two color factors in the quark-annihilation channel which are not yet known analytically. Part of the MIs needed for this task are known from previous works Bonciani:2003te ; Bonciani:2003hc ; Bonciani:2008wf ; Bonciani:2008az ; Bonciani:2009nb ; Bonciani:2010mn ; vonManteuffel:2013uoa ; Bonciani:2013ywa (see also the Loopedia database Bogner:2017xhp ). In particular, the first analytical evaluation of a crossed double box with a massive propagator was presented in vonManteuffel:2013uoa in terms of GPLs. More recently, within the context of a project that requires the analytic evaluation of the NNLO QED corrections to electron-muon scattering Calame:2015fva ; Abbiendi:2016xup , a planar Mastrolia:2017pfy and a crossed DiVita:2018nnh topology, which also enter top-pair production in the $q\bar{q}$ channel, were evaluated analytically using GPLs. In the present work, we provide results for the MIs belonging to the last missing crossed topology and we carry out an independent calculation of the MIs of the topology evaluated in DiVita:2018nnh . These results will allow one to complete the analytic calculation of the two-loop corrections to top-quark pair production in the $q\bar{q}$ channel.

The evaluation of the MIs discussed in this work is carried out by following a by now standard technique based on two steps. First, one observes that the dimensionally regularized scalar integrals which appear in the interference of two-loop and tree-level diagrams can all be written in terms of a reduced set of scalar integrals which are identified as the MIs for the problem under study. The two topologies considered in this work involve 52 and 44 MIs, respectively. The reduction to MIs is carried out by means of the computer programs111Other public programs for the reduction to the MIs can be found in Anastasiou:2004vj ; Lee:2012cn ; Lee:2013mka ; Maierhoefer:2017hyi . FIRE Smirnov:2008iw ; Smirnov:2013dia ; Smirnov:2014hma and Reduze 2 Studerus:2009ye ; vonManteuffel:2012np , that implement integration-by-parts identities Tkachov:1981wb ; Chetyrkin:1981qh ; Laporta:2001dd and Lorentz-invariance identities Gehrmann:1999as . Subsequently, the MIs are computed by employing the differential equations method Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:1999as ; Argeri:2007up ; Henn:2014qga . The system of differential equations is cast in canonical form Henn:2013pwa (see also Argeri:2014qva ; DiVita:2014pza ; Henn:2013nsa ; Gehrmann:2014bfa ; vonManteuffel:2014mva ; Lee:2014ioa ; Adams:2017tga ; Schneider:2016szq ; Meyer:2016slj ; Georgoudis:2016wff ; Gituliar:2017vzm ). The solution is expressed in terms of Chen’s iterated integrals, which can be expanded as a series in the dimensional regularization parameter, and each order of the expansion is represented in terms of GPLs.

The paper is structured as follows. In Section 2, we introduce our notation and we define the topologies that are considered in this work. In Section 3, we briefly review the method of differential equations. In Section 4, we present the canonical form we used for the evaluation of the solution of the system of differential equations. In Section 5, we describe a reparametrization which rationalizes our differential equations. In Section 6, we discuss the integration of the differential equations in terms of GPLs and present the structure of our results. In Section 7, we discuss numerical checks which were carried out in order to validate the analytic expression of the MIs. We emphasize that, in addition to the checks discussed in Section 7, our results have been successfully compared against the expressions of a different set of master integrals, independently obtained by S. Di Vita, T. Gehrmann, S. Laporta, P. Mastrolia, A. Primo, and U. Schubert PD , which were published on the arXiv simultaneously to the present manuscript. Finally, Section 8 contains our conclusions. The definition of the various MIs in terms of momentum integrals over a set of propagators can be found in Appendix A. Numerical results in a specific phase-space point for the seven denominator MIs evaluated analytically in this paper are collected in Appendix B.

Our full analytical results are provided in ancillary files included in the arXiv submission of this paper.

2 Notations

In this paper we consider the process $q\bar{q}\to t\bar{t}$, where $q$ and $\bar{q}$ are massless quarks and $t$ and $\bar{t}$ are massive (top) quarks. The incoming partons have momenta $p_{1}$ and $p_{2}$, while the final state partons have momenta $p_{3}$ and $p_{4}$. All particles are on their mass-shell, namely $p_{1}^{2}\,=\,p_{2}^{2}=0$, and $p_{3}^{2}=p_{4}^{2}=m^{2}$, where $m$ is top-quark mass.

The kinematics of the process can be described in terms of the three Mandelstam invariants

[TABLE]

which satisfy the relation $s+t+u=2m^{2}$. The physical region is defined by

[TABLE]

where $\theta$ is the scattering angle of top quark with respect to the direction of the incoming $q$ quark in the partonic center of mass frame.

Figure 1 shows the two seven-denominator two-loop topologies that we consider in this paper; they are indicated with the capital letters $A$ and $B$. The scalar integrals belonging to Topology $A$ are defined as

[TABLE]

while the scalar integrals belonging to Topology $B$ are defined as

[TABLE]

The labels $a_{i}$ and $b_{i}$, with $i=1,...,9$, are integer numbers where $a_{4}$, $a_{6}$, $b_{5}$, $b_{6}$ $\leq 0$. The $D_{i}$, $i=1,...,9$, are the denominators and numerators involved and $d$ is the dimension of the space-time. The normalization of the integrals is such that

[TABLE]

where $\epsilon=(4-d)/2$, $\gamma_{E}=0.5772..$ is the Euler-Mascheroni constant and $\mu$ is the ’t Hooft scale.

The nine-propagator integral family that we use for the reduction of both Topology $A$ and Topology $B$ is

[TABLE]

Topology $A$ has 52 MIs while Topology $B$ has 44 MIs. Since some MIs are common to both topologies, the total number of independent MIs is 70. Some of the MIs were already known in the literature Bonciani:2003te ; Bonciani:2003hc ; Bonciani:2008wf ; Bonciani:2008az ; Bonciani:2009nb ; Bonciani:2010mn ; vonManteuffel:2013uoa ; Bonciani:2013ywa ; DiVita:2018nnh . Many MIs, including many seven denominator four-point functions, are evaluated here for the first time.

3 The Differential Equations Method

The analytic computation of the MIs is carried out by employing the differential equations method Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:2002zr ; Argeri:2007up ; Henn:2014qga . The MIs can be thought of as components of a vector $\vec{f}$ where each component depends on a vector $\vec{x}$ of dimensionless parameters and on the dimensional regulator $\epsilon$. The dimensionless parameters in $\vec{x}$ are functions of the kinematic invariants of the problem. In the case under study, the vector $\vec{x}$ has two components; the specific choice of these components is discussed in Section 5. The MIs $\vec{f}(\vec{x},\epsilon)$ satisfy a system of first-order partial linear differential equations with respect to the kinematic invariants $\vec{x}$:

[TABLE]

where $\partial_{i}=\partial/\partial x_{i}$, $A_{i}(\vec{x},\epsilon)$ is a set of matrices associated with the system of differential equations. These matrices have dimensions $N\times N$, where $N$ is the number of MIs in the vector $\vec{f}$. In general, the elements of $A_{i}$ also depend on the kinematic invariants, $\vec{x}$, and on the dimensional regulator $\epsilon$. The matrices $A_{i}(\vec{x},\epsilon)$ satisfy the integrability conditions

[TABLE]

where $\left[A_{i},A_{j}\right]\equiv A_{i}A_{j}-A_{j}A_{i}$ is the usual matrix commutator. For a given choice of MIs, the matrices $A_{i}$ can be computed using integration-by-parts identities.

We solve the system in (7) by employing the Canonical Basis approach Henn:2013pwa ; Henn:2014qga , which consists in finding a basis for the MIs in which the system of differential equations has the specific form

[TABLE]

In (9), $d\tilde{A}(\vec{x})$ is a logarithmic differential one-form. Several methods that allow one to find a Canonical Basis for a given topology have been proposed Henn:2013pwa ; Henn:2014qga ; Argeri:2014qva ; Lee:2014ioa ; Georgoudis:2016wff ; Meyer:2017joq . In this work, we find the Canonical Basis by employing the semi-algorithmic approach described in Gehrmann:2014bfa ; Becchetti:2017abb .

In this basis the solution of the system of differential equations in (9) is formally written in terms of Chen iterated integrals Chen:1977oja :

[TABLE]

where $\mathbb{P}$ stands for the path-ordered integration, $\gamma$ is some path in the space of kinematic invariants and $\vec{f}_{0}(\epsilon)$ is a vector of boundary conditions that we found by imposing the regularity of the MIs in particular points of the phase space or known solutions for simple integrals.

For the process under study it is possible to find a change of variables such that the matrix $\tilde{A}(\vec{x})$ is a rational function of the kinematic invariants, i.e. the entries of the one-form $d\tilde{A}(\vec{x})$ are linear combinations of terms $d\log(x_{k}-\alpha_{k})$, where $\alpha_{k}$ are algebraic functions of kinematic invariants and the arguments of the logarithms $(x_{k}-\alpha_{k})$ determine the so called alphabet of the process. In this case, once a path $\gamma=\gamma(t)$ is specified, the solution can be written order-by-order in the dimensional regularization parameter explicitly in terms of GPLs

[TABLE]

with

[TABLE]

where $\vec{0}_{n}$ indicates a list of $n$ weights, all equal to [math].

4 Canonical Form for the Master Integrals

In this Section, we present the canonical basis used in this work. In particular, we provide the relations that allow one to go from the MIs in pre-canonical form (see Fig. 2) to MIs in canonical form, where the latter satisfy a differential equation of the form (9). These relations are written assuming that the Mandelstam invariants $s$ and $t$ take values in the physical region of phase space. The normalizing prefactors contain two square roots, one of which enters only through the two integrals $f^{A}_{44}$ and $f^{B}_{38}$, see below. As will be shown in the next section, it is possible to rationalize these roots by a suitable reparametrization. We will extend the definition of the canonical MIs also to other phase-space regions by using the expressions listed in (13 – 108) after rationalization. In other words, we effectively define our canonical MIs using rational prefactors in the parameters {$w$, $z$} rather than analytically continue root-valued prefactors in the original Mandelstam invariants.

The pre-canonical MIs basis is shown graphically in Fig. 2. As discussed in the caption, thin lines represent massless propagators and external legs, while massive lines represent massive propagators and external legs. For two and three point functions, we indicate explicitly the Mandelstam variable on which a given MI depends by adding “$s$”, “$t$”, “$u$” or “$p_{3}^{2}$” to the drawing (see for example ${\mathcal{T}}^{A}_{9}$,${\mathcal{T}}^{B}_{2}$, etc.). For sub-topologies involving several MIs, dotted propagators indicate a squared propagator in the integrand of the MI (see for example ${\mathcal{T}}^{B}_{17}$). A dotted propagator with a $3$ next to the dot indicates a cubed propagator in the integrand (see for example ${\mathcal{T}}^{B}_{12}$). The four-point subtopologies in the last two lines of Fig. 2 involve several MIs which differ from one another because of the numerator in the integrand. These numerators are shown on top of the drawing of each single MI (see for example ${\mathcal{T}}^{B}_{37}$, etc). In order to avoid any possible misinterpretation of Fig. 2, the integral definition of each MIs in the pre-canonical basis can be found in Appendix A. In addition, we provide all definitions in the ancillary files of the arXiv submission of this work.

4.1 Topology $A$

Topology $A$ involves 52 MIs. Their canonical form is obtained with the following change of basis:

[TABLE]

4.2 Topology $B$

Topology $B$ involves 44 MIs. The relations linking the canonical and pre-canonical forms of the MI basis are the following:

[TABLE]

5 Rationalizing Parametrization

Both the definition of the Canonical Basis for topologies $A$ and $B$ and the differential equations fulfilled by the MIs in the Canonical Basis contain two square roots of functions of Mandelstam invariants. These roots are:

[TABLE]

where the latter root enters only through the definition of two MIs in the canonical basis: $f^{A}_{44}$ and $f^{B}_{38}$. In order to express the solution in terms of GPLs one needs to rationalize the square roots by finding an appropriate change of variables. We start by defining the dimensionless variables $x$ and $y$ as follows:

[TABLE]

A particularly convenient parametrization LorenzoReparam:2018 ; DiVita:2018nnh which rationalizes the roots222Another solution to the rationalization problem can be obtained from diophantine equations as described in Becchetti:2017abb : $x={16w^{2}(1+4z)^{2}(w+z+4wz)^{2}}/[{z(1+8w)(z^{2}-4w^{2})(z+4w+8wz)}]$ , $y=8z(1+2z)$. However, this change of variables leads to rather long expressions. in (109) is given by

[TABLE]

Since we want to express all of the MIs in terms of $w$ and $z$, it is also necessary to write the Mandelstam variable $u$ in terms of $w$ and $z$. Starting from the relation $u=-s-t+2m^{2}$ one finds

[TABLE]

The analytic expressions that can be found in the ancillary files are written in terms of $w$ and $z$ defined in (111).

In order to write $w$ and $z$ as a function of $x$ and $y$, it is necessary to invert the system in (111). The first equation has two solutions

[TABLE]

When $0<x<\infty$ ($-\infty<s<0$), the first solution is limited and such that $0<w_{1}<1$, while the second solution is unlimited, $1<w_{2}<\infty$. In view of the fact that GPLs of $w$ should be manifestly real in this region, we choose the first solution, $w_{1}$. The second equation, (111), has four solutions, two for each choice of $w$:

[TABLE]

The first solution $z_{1}$ is always positive for $y>0$ ($t<0$) and $w=w_{1}$, while the second is always negative. When, for a given $x$, $y$ varies from 0 to $\infty$, $z_{1}$ is limited to the range $\sqrt{w}<z_{1}<\sqrt{1-w+w^{2}}$.

Consequently, in the region $s<0$ and $t<0$, we choose the set of variables

[TABLE]

When $s$ becomes positive, it is necessary to consider the Feynman prescription and add to $s$ a positive vanishing imaginary part: $s+i0^{+}$:

[TABLE]

where, now, $x^{\prime}=s/m^{2}>0$. If $x^{\prime}$ is such that $0<x^{\prime}<4$, $w$ becomes a phase and it moves on the upper unit circle:

[TABLE]

with

[TABLE]

For $x^{\prime}=4$, $w$ becomes real again and one finds that $w=-1$.

For physical kinematics one finds that $s>4m^{2}$, $t<0$, $u<0$, where

[TABLE]

with

[TABLE]

In this physical region we use

[TABLE]

with the phase space constraint

[TABLE]

The crossing $t\leftrightarrow u$ is given by

[TABLE]

By keeping into account the relation (123), the roots in (109) become

[TABLE]

We apply these substitutions to the definition of the canonical integrals in the physical region in equations (13 – 108) and employ the resulting rational expressions to define the Canonical Basis also in other regions of the phase space.

6 Integration and Results

In terms of the rationalizing variables, we can write for the matrix $\tilde{A}(\vec{x})$ in our differential equation (9)

[TABLE]

where the $\tilde{A}^{(k)}$ are rational matrices and the letters $l_{k}$ form the alphabet

[TABLE]

The last letter is needed only for Topology B. We provide explicit expressions for the matrix $\tilde{A}(\vec{x})$ for Topologies $A$ and $B$, respectively, in the ancillary files on arXiv.

This alphabet allows to analytically integrate the MIs in terms of GPLs of argument $w$ with the weights

[TABLE]

and GPLs of argument $z$ with the weights

[TABLE]

We fix the boundary constants by imposing regularity conditions, supplemented by external input for a few well-known simple integrals.

The analytic continuation of the GPL functions of $w$ and $z$ between different regions of the phase space is non-trivial. We found it convenient to provide the MIs in terms of an analytic expression which is valid in the region $s<0$ and of a second analytic expression that is valid in the region $s>0$. Our complete results are available in the ancillary files sol-A-unphys.m, sol-B-unphys.m for $s<0$ and sol-A-phys.m, sol-B-phys.m for $s>0$.

7 Numerical Checks

In order to validate our results we performed numerical checks in different points of the phase space. In several cases, the MIs were checked by evaluating the MIs in the pre-canonical basis numerically by means of Sector Decomposition Binoth:2000ps as implemented in SecDec Carter:2010hi ; Borowka:2012yc ; Borowka:2015mxa ; Borowka:2017idc and FIESTA Smirnov:2008py ; Smirnov:2013eza ; Smirnov:2015mct and by subsequently comparing the numerical results with the evluation of the analytic expressions for the MIs carried out with GiNaC Bauer:2000cp .

However, in some cases, and in particular for the MIs involving six or seven denominators, we were not able to obtain sufficiently precise numbers by a direct evaluation of the MIs with SecDec or FIESTA. For this reason, we employed the techniques described in vonManteuffel:2014qoa ; Panzer:2014gra ; vonManteuffel:2017myy and rewrote the canonical MIs as linear combinations of quasi-finite integrals. Quasi-finite integrals are integrals which have, at worst, a single pole in $\epsilon$ which originates from the Euler Gamma function prefactor in the Feynman parameter representation of the integral. Quasi-finite integrals are built with the same set of propagators as the original integral but they might be defined in shifted space-time dimensions and might have one or more propagators squared or raised to higher power. It was shown in vonManteuffel:2014qoa ; vonManteuffel:2017myy that quasi-finite integrals are evaluated more efficiently by SecDec and FIESTA with respect to non quasi-finite integrals with the same sets of propagators. Using SecDec we generated numerical results for quasi-finite integrals in the unphysical and in the physical region. Subsequently, we converted these numbers to results for the canonical integrals, at which level we were left with typically 2-6 significant digits, depending on the integral and the region of phase space. We successfully compared these numbers against those obtained from the analytic expressions of the MIs, which are the main result of the present work. With this procedure it was possible to test numerically all of the 52+44 MIs evaluated in this paper.

In addition, we compared numerically the MIs that are in common with the ones presented in DiVita:2018nnh with a numeric evaluation of their expressions, finding complete agreement. Finally, all of the MIs evaluated in this work were simultaneously evaluated in PD . We compared numerically the MIs evaluated in this work with the results obtained in PD . This comparison was carried out in several phase space points, both in the physical region ($s>4m^{2}$) and in the non-physical region ($s<0$). We found complete agreement between the results in this work and the ones in PD .

8 Conclusions

In this paper we presented the analytic calculation of the master integrals necessary for the evaluation of the last two color coefficients of the interference between two-loop and tree-level diagrams for the partonic process $q\bar{q}\to t\bar{t}$, for which an analytic expression is not yet available.

The master integrals were evaluated with the method of differential equations. By determining a canonical basis, we brought the system of first-order linear differential equations into an $\epsilon$-form, allowing for their decoupling after an expansion in powers of the dimensional regulator $\epsilon$. We integrated the expansion coefficients in terms of generalized harmonic polylogarithms of two dimensionless variables through to weight 4 and fixed the integration constants using regularity conditions and known solutions for simple integrals.

We checked our analytic results numerically against the results obtained with two numerical codes, SecDec and FIESTA, using the method of quasi-finite integrals. We also compared numerically the MIs that are in common with the ones presented in DiVita:2018nnh , finding complete agreement. Finally, we cross-checked numerically the MIs of our topology A with the authors of PD , in several points of the phase space, finding complete agreement.

All the analytic expressions for the MIs presented in this paper are provided as computer readable ancillary files together with the arXiv submission of this work.

9 Acknowledgments

We are grateful to L. Tancredi for early collaborations and for providing the rationalizing parametrization used in this work. We would like to thank S. Di Vita, T. Gehrmann, S. Laporta, P. Mastrolia, A. Primo, and U. Schubert for the cross checks of the MIs presented in this work, in several points of the phase space. The work of A.F. is supported in part by the National Science Foundation under Grant No. PHY-1417354 and PSC CUNY Research Award TRADA-61151-00 49. The work of A.v.M. is supported in part by the National Science Foundation under Grant No. 1719863. Feynman diagrams are drawn with Axodraw Vermaseren:1994je .

Appendix A Pre-canonical form for the Master Integrals

In this section we present the routing for all the pre-canonical MIs of the two topologies. The definition of the integration measure ${\mathcal{D}}^{d}k_{i}$ ($i=1,2$) can be found in (5). The list of the denominators $D_{i}$ ($i=1,\cdots,9$) can be found in (6).

A.1 Topology $A$

[TABLE]

A.2 Topology $B$

[TABLE]

Appendix B Numerical Results

In this Appendix we collect numerical results for the seven-denominator canonical MIs at the point

[TABLE]

The numerical values of the MIs are (with 16 significant digits):

[TABLE]