Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel

We present a complete set of analytic helicity amplitudes for top quark pair production via gluon fusion at two-loops in QCD. For the first time, we include corrections due to massive fermion loops which give rise to integrals over elliptic curves. We present the results of the missing master integrals needed to compute the amplitude and obtain an analytic form for the finite remainders in terms of iterated integrals using rationalised kinematics and finite field sampling. We also study the numerical evaluation of the iterated integrals.


Introduction
Precise understanding of the top quark in the Standard Model is one of the highest priorities in modern collider experiments. The large mass of the top quark makes it a dominant ingredient in understanding the fundamental forces and in particular of the electroweak symmetry breaking mechanism. A precise determination of the top quark mass is imperative due to the sensitivity of the Higgs potential on its value.
The QCD corrections to top-quark pair production have a long history going back to the first total cross section computation at NLO for on-shell top-quarks computed by Nason, Dawson and Ellis [1,2]. Thanks to a monumental effort on the part of the theoretical community, fully differential NNLO predictions for top anti-top pair production at hadron colliders are now available for comparisons with the latest experimental data [3][4][5][6][7]. These intensive computations require two-loop 2 → 2, one-loop 2 → 3 and tree-level 2 → 4 partonic amplitudes to be combined using a consistent UV and IR subtraction scheme to obtain a finite result. The first complete computation used the 'STRIPPER' sector decomposition method to subtract IR divergences and has enabled the most complete theoretical descriptions of top-quark dynamics [8]. Very recently a second independent computation has been obtained using the q T subtraction method [9][10][11][12].
Up to now the only complete two-loop amplitudes have been available numerically [13][14][15]. The main reason for this is that a more complicated class of special functions begin to appear in amplitudes which contain internal masses. These functions have been identified to involve integrals over elliptic curves [16][17][18][19][20][21][22][23][24], and lie on the boundary of our current mathematical understanding. Over the last few years, motivated by the increasing demand of precise theoretical descriptions including mass dependence, there has been substantial progress in developing a complete framework for use in phenomenological applications. There are a large number of amplitude level corrections to the process which do not depend on elliptic sectors and there has been substantial effort to find a complete analytic form for the squared amplitudes [25][26][27][28][29][30][31][32].
In this paper we present a set of helicity amplitudes for top-quark pair production in the leading colour approximation. For the first time we also include contributions from heavy fermion loops which lead to the appearance of iterated integrals involving elliptic curves. The master integrals for these heavy loop corrections have been recently completed [17,18]. As well as obtaining compact analytic helicity amplitudes by sampling Feynman diagrams with finite field arithmetic [33][34][35][36][37][38], we also study the numerical evaluation of the final amplitudes. The helicity amplitudes presented contain complete information about top quark decays in the narrow width approximation. This helicity amplitude technique has been used successfully at one-loop [39][40][41][42]. Spin correlations can also be treated efficiently using the spin density matrix approach [15,43,44] As is always the case, computations at this perturbative order contain a large number of steps, each usually with some technical bottlenecks to overcome. The analytic computations of the planar master integrals, including those with an internal massive loop, are available in the literature for almost every topology. Yet, for the construction of the amplitude it was necessary to include one final two-loop integral topology which contained two elliptic master integrals. The (canonical form) differential equation [21,[45][46][47][48][49] and solution in terms of iterated integrals was obtained and combined with the integration-by-parts reduced helicity amplitudes.
Throughout the computation we made use of finite field arithmetic to find an efficient solution to the system of integration-by-parts identities. These techniques have shown to be particularly efficient for massless amplitudes with many external scales. In this paper we show how they can apply equally well to amplitudes with massive internal particles. We construct a rational parametrisation of the on-shell kinematics via momentum twistors and obtain helicity amplitudes via projecting onto a basis of independent spinor structures, accounting for the freedom in the top quark spin states. This method defines a set of on-shell, gauge invariant sub-amplitudes that can be computed using on-shell top-quark kinematics. We then check these results against previous computations for the interference with tree-level diagrams in the conventional dimensional regularisation (CDR) scheme.

Leading colourttgg amplitudes
We consider a scattering process involving a pair of top quarks and two gluons 0 →t(p 1 ) + t(p 2 ) + g(p 3 ) + g(p 4 ), where p 2 1 = p 2 2 = m 2 t and p 2 3 = p 3 4 = 0. The kinematic invariants for this process are the top-quark mass m t , and the two Mandelstam variables In this work we consider the leading colour contributions of thettgg amplitude up to two-loop level, where at two loops, only planar configurations arise. The colour decomposition of the leading colour L-loopttgg amplitude is given by where n = m α s /(4π), α s = g 2 s /(4π), m = i 4π/m 2 t e − γ E , g s is the strong coupling constant and (T a )j i are the fundamental generators of SU (N c ). The partial amplitudes can further be decomposed according to their internal flavour structure, where N l and N h are the number of closed light and heavy quark loops, respectively. Sample diagrams for various fermion loop contributions at one and two loops are shown in Figs. 1 and 2. We first define the gauge invariant L-loop mass-renormalised amplitude A (L) mren , where the top-quark mass renormalisation is included mct ) represents an amplitude including one mass counterterm insertion to top-quark propagators in the tree level (one-loop) diagrams, while A (0) 2mct corresponds to the tree level amplitude with two mass counterterm insertions to top-quark propagators. The mass counterterm insertion procedure leads to the presence of double and triple top-quark propagators in the tree level amplitude and double top-quark propagators in the one-loop amplitude. Sample diagrams with mass counterterm insertions at tree level and one loop are displayed in Fig. 3.
The fully renormalisedttgg amplitude can be obtained by including renormalisations of the top-quark and gluon wavefunctions as well as the strong coupling α s and δZ (L) αs are the L-loop top-quark wavefunction, gluon wavefunction and strong coupling renormalisation constants, respectively. All the renormalisation constants relevant for our calculation are given in Appendix A.
The renormalisation procedure removes the ultraviolet (UV) singularities present in loop amplitudes. The remaining divergences that are of infrared (IR) nature can be predicted from universal IR behaviour of QCD amplitudes [50][51][52][53][54], which, in the presence of top quarks, needs to be extended to massive case [55][56][57]. The UV renormalised amplitude at one and two loops can be divided into singular and finite terms (2.9) (2.10) Here |A (L) n is a vector in colour space and admits the following expansion in the strong coupling We remind that for the leading colourttgg amplitudes the colour decomposition is given in Eq. (2.2). The IR divergences of the amplitude are encoded in the Z factor (which is a matrix in colour space) and up to two-loop order. For the process we are considering, the Z factor is given by [57] In this work we will consider the top quark as the only massive fermion, therefore n h = 1 and m 1 = m t . The coefficients Γ n and Γ n are defined through the following expansion in α s , For thettgg amplitude, C i = C A = N c . The anomalous dimension matrix Γ for thettgg amplitude can be obtained from Eq. (55) of [57]. We also refer the reader to Appendix A of [52] for the definition of anomalous dimension factors required for the computation of the IR poles of thettgg amplitude.

Helicity amplitudes for massive fermions
There are a variety of different methods on the market for dealing with the spin states of a massive fermion. The common aim is to find a compact notation at the amplitude level and retain all information required to account for decays in the narrow width approximation.
One major issue is that the helicity for a massive particle is not conserved and so we must introduce an arbitrary additional direction in order to define it. This additional extra direction can increase the algebraic complexity of an analytic computation.
In this section we review one way in which massive spinors can be incorporated into the spinor-helicity methods [58] and then describe how to define a set of gauge invariant, onshell sub-amplitudes that can be used to describe the full set of spin correlated narrow width decays. These on-shell sub-amplitudes can be computed using a rational parametrisation of the external kinematics. In our case, this rational paramtrisation was generated using the momentum twistor formalism [59]. The spinor-helicity method applied to massive fermions has been well studied and we refer the reader to Refs. [60][61][62] for further details and other approaches.
The first step is to introduce an arbitrary direction n which can be used to define a massless projection of the massive fermion momentum, where p µ is the momentum of an on-shell massive fermion with p 2 = m 2 . The momenta p and n are both massless. The set of u and v spinors for the massive fermion can then be constructed using the Weyl spinors for the massless momenta, The fact that the helicity state depends on the choice of reference vector means that positive and negative helicities are no longer independent as they are in the massless case: For a tt pair this means that we only need to compute one helicity configuration with general reference vectors, say ++, and we will obtain enough information for all four possible spin configurations.
We can then organise amplitudes for the ++ configuration into a basis of spinor structures which parametrise the dependence on the two reference vectors n 1 and n 2 . There are four independent terms in this basis, matching up with the four configurations in the complete system. For the L-loopttgg amplitude this is where Φ is a phase factor depending on the helicities of the gluons. There are a few things to note about this expression: 1. The sub-amplitudes A (L), [a] are defined to be dimensionless and free of the spinor phases.

A (L),[a]
depends only on two variables.
3. The gluon phases for the ++ and +− configuration are taken to be 4. The basis of spin structures has been chosen to prefer n 1 n 2 over n 1 4 n 2 3 . The sub-amplitudes appear to be simpler in this basis.
The computation of the sub-amplitudes can then be obtained from four different evaluations of the full amplitude with a rational kinematic configuration with four choices of the reference vectors. These evaluations can then be used to make a linear system which is solved to find the sub-amplitudes A (L), [a] . Explicitly, A (L), [4] , (3.6) A (L), [3] , (3.7)

8)
where we have dropped the particle labels on the A functions for simplicity.

Generating a rational parametrisation of the kinematics
We start from a rational configuration with six massless particles generated via momentum twistors [59,63]. This configuration will depend on eight independent parameters. The six particles are then used to form a configuration with specific choices for the reference momenta. We label the six massless momenta as q 1 , . . . , q 6 and the resulting tt system as p 1 , p 2 , p 3 , p 4 where, The ordering of the momenta q 1 , . . . , q 6 is important to find a simple rational solution to the last four constraints. From the on-shell momenta p 1 , . . . , p 4 it is straightforward to generate rational spinors for p 1 , p 2 using four different choices of reference vectors: (n 1 , n 2 ) = {(3, 3), (4,4), (3,4), (4, 3)}. The variables of the original six-point configuration, i.e. those of the q i momenta, are changed so that we can use conventional Mandelstam invariants, This procedure would also work for high multiplicity amplitudes though a careful choice of variables in the parametrisation may be necessary to obtain manageable algebraic complexity.
The final results for the sub-amplitudes in Eq. (3.5) will be expressed using kinematic variables x and y defined by This choice rationalises the square root, that appears in the master integral and amplitude computations. For thettgg amplitudes that do not involve elliptic sectors, the rationalisation of such a square root allows us to perform a computation within finite field arithmetic without the need to introduce additional variables. It also important to make a rational change of variables when solving the differential equations of the master integrals.

Tree-level sub-amplitudes
Following the procedure above allows us to directly evaluate the relevant colour ordered Feynman diagrams to obtain simple expressions at tree level: [2] (1 + t , 2 + The procedure for generating loop level expressions does not change although the expressions require, as usual, tensor integral reduction to be performed. Our implementation of these steps is given in the subject of the following sections. Our final aim is to present the finite remainders of the sub-amplitudes after the subtraction of infrared and ultraviolet poles.

Amplitude reduction
The reduction of the helicity amplitudes, generated using Feynman diagrams, is performed directly to special functions using finite field arithmetic in the FiniteFlow framework [38]. Recent years have seen a growth in the popularity of loop amplitude computations with finite field arithmetic for problems with many external scales where the algebraic complexity is too high for conventional approaches. There are however many potential bottlenecks in amplitude reduction and it has been shown to be important to control as many steps in the computation as possible. Amplitudes requiring dimensional regularisation can cause unnecessary large intermediate steps since information that vanishes in four dimensions is retained. Unitarity cut based approaches [64] and specially designed tensor projectors [65][66][67] have been used successfully in recent high multiplicity loop amplitude computations in massless theories. In our setup the diagram generation, colour ordering and spinor-helicity algebra are performed with the help of Qgraf [68], Form [69,70], Spinney [71] and Mathematica scripts. During this phase we identify a set of topologies that are independent at the integrand level. The numerators of these topologies are then separated into loop momentum dependent structures and coefficients depending on the external kinematics. The latter are then evaluated using the rational parametrisation described in the previous section. Up to this point we follow a strategy described in, for example, Ref. [72]. The only difference is that the internal masses must also be tracked.
The result of these steps is not yet suitable for processing with integration-by-parts identities since the loop dependent numerators must be transformed into a basis of independent scalar products thereby upgrading the topologies to complete families of Feynman integrals. In the case of four-particle amplitudes an additional integration over the spurious space in the loop momenta must also be performed. We employ a transverse integration approach similar to the one outlined in Ref. [73].
For clarity, we present some additional details of these steps that are specific to our problem. An interested reader may wish to refer to Refs. [74][75][76] for a more complete discussion and other applications. All steps described from this point to the reduction onto a basis of special functions have been constructed using FiniteFlow graphs. After colour ordering and helicity amplitude processing the amplitude takes the following form, where d = 4 − 2 , d s = g µ µ and k i are the spacetime dimension, spin dimension and loop momenta, respectively. T is a set of distinct diagram topologies contributing to the colour- Since the external momenta are in four dimensions, the extra dimensional part of the loop momenta only appears in the numerator function as µ ij = −k i ·k j . After substituting massive spinors in terms of massless ones according to Eqs. (3.2) -(3.3) and performing t'Hooft algebra, the numerator function contains the following loop-momentum dependent terms where momenta p a and p b are massless and p j can either be massive or massless. As indicated above, to apply IBP reduction we must first express these loop momentum structures in a basis of the propagators D α and a set of irreducible scalar products (ISPs) which define a family of Feynman integrals. We make sure that we introduce the minimal number of integral families by mapping all distinct diagram topologies, T , to a set of master topologies which have the maximum number of propagators. The master topologies for the N 2 c part of the two-loop planarttgg amplitude are shown in Fig. 4. Each topology T will map to at least one master topology. The numerators are then expressed in terms of the loop momenta of the master topology. We split the four-dimensional loop momenta k i into physical partsk µ ,i and spurious parts and we further expand them into a physical (v µ ) and a spurious (w µ ) spanning bases is the dimension of the 4d physical space, n ext is the number of independent momenta of a given topology and d ⊥ = 4 − d is the dimension of the 4d spurious space. In addition, the physical vectors are orthogonal to the spurious vectors, v i · w j = 0.
The coefficients of the physical spanning basis are functions of inverse propagators and physical ISPs, a ij = a ij (D α , ISPs), while the coefficients of the spurious spanning basis are functions of spurious ISPs and extra dimensional ISPs . a ij and b ij are determined by solving a linear system of equations obtained by contracting Eq. (4.5) with v i and w i . For thettgg amplitudes, all the maximal topologies have four external momenta, hence n ext = 3, d = 3 and d ⊥ = 1, and for instance, we can choose the following set of Explicitly, we used: though any vector satisfying ω · p 1 = ω · p 3 = ω · p 4 = 0 will suffice. Numerator terms with odd powers of k i · w j vanish under loop integration, while the terms with even powers can be cast into linear combination of propagator denominators D α and physical ISPs. For example, where k i · k j can be directly written in terms of D α while k ,i · k ,j can be expressed in terms of D α and physical ISPs according to Eq. (4.5). Similarly for the extra dimensional ISP, After the replacements have been made the amplitude is recast into a form suitable for IBP reduction, Note that there may be overlap in sets of Feynman integrals G T,i ({k}, {p}) appearing in each numerator topology T . This is in contrast to the integrand reduction approach taken in [72] where there is no overlap in G T,i ({k}, {p}) between numerator topologies. Combining contributions from all numerator topologies we have whereG k is a unique set of Feynman integrals obtained from all the G T,k in Eq. (4.10).
Before proceeding with IBP reduction, at this stage we combine the bare helicity amplitude with the mass renormalisation counterterm contributions according to Eqs. (2.5) and (2.6), for the one-and two-loop amplitudes, to obtain a gauge invariant, mass renormalised amplitude, The mass renormalisation counterterms appearing in Eqs. (2.5) and (2.6) have to be written in the integral representation, in order to be combined with the bare amplitude. The integrals appearing in Eq. (4.13) are further reduced to a set of master integrals via IBP identities. The IBP relations are generated in Mathematica with the help of LiteRed [77], and solved numerically over finite fields within the FiniteFlow framework [38] using the Laporta algorithm [78]. The mass renormalised helicity amplitude is now written in a basis of master integrals MI k , (4.14) Let us note that not all of the maximal topologies that appear in an amplitude are processed in the IBP system, as some of them can be mapped into the other maximal topologies. For example, maximal topologies T 1 1 , . . . , T 1 4 in Fig. 4 can be mapped to topologies T 1 5 , . . . , T 1 8 by choosing an appropriate set of physical ISPs (auxiliary propagators) for each of those topologies. The master integral representation can be used to check Ward identities on the mass-renormalised amplitude by replacing gluon polarisation with its momentum in the numerator construction and reconstructing Eq. (4.14) over a prime field to check if it indeed vanishes.
If the analytic solutions of the master integrals appearing in Eq. (4.14) for the planarttgg amplitudes are available, we can proceed further by writing the master integrals as a linear combination of special function monomials m (f ({p})) and perform Laurent expansions, where n(L) is the power of the deepest pole that can appear in the L-loop amplitude. For one-and two-loop cases, n(1) = −2 and n(2) = −4. To obtain the simplest representation of the amplitude we further subtract the UV and IR poles from the mass-renormalised amplitude. We first need to write the UV and IR poles, given in Eqs. (2.7), (2.8), (2.9) and (2.10) for the case of planarttgg one-and two-loop amplitudes, in the same special function monomial basis as in Eq. (4.15), and subtract them from the mass-renormalised amplitude to obtain the finite remainder, The coefficients c F,h k ({p}) appearing in (4.19), however, are not all independent. We exploit this fact in order to simplify both the result and the reconstruction of its analytic expression as follows. First, we sort all the coefficients by their complexity, which is estimated from their total degree. The total degree, in turn, can be quickly determined via a univariate fit, as explained in [37]. We then find vanishing linear combinations of these coefficients by solving the linear fit problem with respect to the unknowns y k , which are rational numbers. This allows us to find linear relations between the coefficients, which rewrite the more complex ones in terms of simpler ones. After applying the linear relations between coefficients of special function monomials, we arrive at wherec F,h k are the independent coefficients of the new special function monomialsm k (f ), which are linear combinations of the monomials appearing in Eq. (4.19). Therefore, functional reconstruction only needs to be applied to the independent coefficientsc F,h k . This yields a significantly simpler result than the one in Eq. (4.19), and also reduces the number of evaluations needed for its reconstruction.
Having set up the numerical algorithm suitable for computation over finite fields, starting from the Feynman diagram numerators, to evaluatec F,h k in the L-loop finite remainder F (L),h n ({p}) using the FiniteFlow framework, the analytic forms ofc F,h k are obtained by using the multivariate reconstruction method described in [37], after performing several numerical evaluations.

Master integrals
To derive an analytic representation of the helicity amplitudes, particularly when we are interested in computing the finite remainder, analytic expressions of master integrals appearing in the one-and two-loop amplitudes are required. The solutions of all one-and two-loop master integrals appearing in the A (2),1 , A (2),N l , A (2),N 2 l , A (2),N l N h and A (2),N 2 h amplitudes in Eq. (2.4) can be expressed in terms of multiple polylogarithms (MPLs) [79][80][81] and have been completed recently [27,82]. For the two-loop amplitude A (2),N h involving a single top-quark closed-loop, the master integrals also involve elliptic generalisations of MPLs [17,18]. Before we proceed further, it is useful to introduce the notion of a sector (or a topology). A sector is defined by the set of propagators with positive exponents. In our work, we employ the following master integrals to obtain analytic helicity amplitudes:  • A (2),1 , A (2),N l , A (2),N 2 l : There are 36 master integrals appearing in the A (2),1 amplitude. We have used the results in [82] for the relevant master integrals. The basis of master integrals involves the topologies with 7 propagators (double-box) shown in Fig. 5 (topology a and b). The master integrals required for the computation of A (2),N l and A (2),N 2 l can be obtained from a subset of the master integral basis of A (2),1 .
• A (2),N h : There are 55 master integrals in the amplitude with a single top-quark closed loop. 44 of them belong to the topbox family, i.e. the planar double-box integral family contributing to A (2),N h (topology c in Fig. 5), computed in [18]. Results for the 9 master integrals that are not part of the topbox family are obtained from [82]. Analytic expressions for the remaining 2 master integrals (shown in Fig. 6), to the best of our knowledge, are not available in the literature. In this work we have computed those two integrals using the method of differential equations. The details of the computation will be covered in the next subsection.
• A (2),N l N h and A (2),N 2 h : Both amplitudes are of one-loop squared type, and the master integrals needed are the one-loop bubbles and triangle up to weight 4, which are already available from the one-loop amplitudes computation up to O( 2 ). Figure 7: Feynman diagram contributing to A (2),N h that leads to the new master integral topologies shown in Fig. 6. The solid lines denoted the top quark whereas spiral lines denote the gluons. All external momenta are on-shell.

New master integrals
The new master integral topologies that appear in A (2),N h arise from the Feynman diagram shown in Fig. 7. The integral family for this penta-triangle topology is given by where γ E denotes the Euler-Mascheroni constant, ν = 9 j=1 ν j and d = 4 − 2 . The inverse propagators D i are given by We choose a basis of master integrals for which the irreducible scalar products corresponding to D 8 and D 9 are absent. Therefore we may label these master integrals by I ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 .
Using IBP identities, we obtain 22 master integrals, {M i } 22 i=1 which are shown in Fig. 8. We refer to this basis of master integrals as the pre-canonical basis, and they are given by where I ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 are the precanonical master integrals defined in Eq. (8) of [18]. A second set of 5 master integrals is available from [82], and they are identified as where T i are master integrals defined in Sec. 6.2 of [82]. M 21 and M 22 are the missing integrals. Note that these two integrals also appear in a planar double-box family that contributes to the subleading colourttgg amplitude without closed fermion loops (topology P 5 in Figure 1 of [83]). Apart from the x and y defined in Eq. (3.13), it can be useful to use another coordinate systemx andỹ defined as in order to rationalize some of the square roots appearing in the differential equations.
In the following we aim to bring all the master integrals in a canonical form [21,49]. For the master integrals not belonging to the topbox family we change the basis such that where A does not depend on and is rational in the kinematic variables (x, y). In particular, the canonical basis for integrals in Eq. (5.5) along with the missing integrals from that family is chosen as For the set of master integrals belonging to the topbox topology we follow the choice given in Sec. 6 of [18]. Here the form of the differential equation is slightly relaxed and the differential equation for the basis J has the form where A (0) is strictly lower triangular and A (0) and A (1) are independent of . A (i) x,y is rational in the kinematic variables (x, y), periods of elliptic curves and their derivatives [21]. By considering them as independent variables, we derived the system of differential equations efficiently with the finite field sampling method using the FiniteFlow package [38]. We note that in the computation of the A (2),N h finite remainder, the coefficients of the special function monomials,c F,h k in Eq. (4.21), are also rational functions in (x, y), periods of elliptic curves and their derivatives. The solution for master integrals in terms of iterated integrals are presented up to 4 orders in : (5.10)

Iterated integrals
Now that we have identified the necessary master integrals, let us discuss the mathematical objects that can used to describe them. By iteratively solving the differential equation (5.9), we naturally obtain a representation of master integrals in terms of iterated integrals. Let us therefore start by reviewing the definition of Chen's iterated integrals [84]. Consider a path γ on an n-dimensional manifold M with starting point x i = γ(0) and end point Let us also consider a set of differential 1-forms {ω i } and their pullbacks to the interval [0, 1] which we denote by Then we define the k-fold iterated integral over ω 1 , ..., ω k along γ as From the above definition it is clear that the derivative of an iterated integral is given by d dλ I γ (ω 1 , ω 2 , ..., ω k ; λ) = f 1 (λ) I γ (ω 2 , ..., ω k ; λ) .
It can be shown that the iterated integrals defined in Eq. (5.12) have various interesting properties [85]. The most useful property for the purpose of this paper is concerning the decomposition of the path γ. Let α, β: [0, 1] → M be two paths such that β(0) = α(1) and let γ = αβ be the path obtained by concatenating α and β. Then we can decompose I γ (ω 1 , ..., ω k ; λ) = n i=0 I β (ω 1 , ..., ω i ; λ) I α (ω i+1 , ..., ω n ; λ) , (5.14) where we define the 0-fold integrals as Now that we are familiar with the notion of iterated integrals, let us consider the wellunderstood multiple polylogarithms (MPLs). MPLs are a special class of iterated integrals where the 1-forms ω i are such that for some c ∈ C. Then we have Note that the recursion defined in (5.17) diverges whenever z k = 0. We can regularise this divergence by setting The class of MPLs is very well understood and many tools for their algebraic manipulation and their numerical evaluation are available [86][87][88], it is however not always possible to express all master integrals in terms of these functions. Especially in computations involving massive internal particles we often encounter elliptic sectors which transcend the space of MPLs. In the following, we give a brief overview of the elliptic sectors we encounter in this work as well as the computational peculiarities that come along with them.

Elliptic curves
In this section we review some basic properties of elliptic curves as well as the elliptic curves specific to the topbox topology, which were previously analysed in [18]. For more thorough reviews of elliptic curves and their properties, we refer the reader to one of the various analyses in the literature [18,89,90]. Let us start by considering the generic quartic case of an elliptic curve. An elliptic curve E can be described by the equation where the z i ∈ C. The z i define the properties of the elliptic curve and are generally functions of the kinematic variables x = (x 1 , ..., x n ): Alternatively, an elliptic curve can be described in terms of two elliptic periods ψ i . In order to define these elliptic periods, let us introduce the auxiliary variables Z i as The Z i satisfy the relation and can be used to define the modulus k 2 and the complementary modulusk 2 of the elliptic curve E as Then we can choose the two elliptic periods ψ i associated to the elliptic curve E as where K denotes the complete elliptic integral of the first kind In our results, the dependence of iterated integrals on elliptic curves enters through the appearance of elliptic periods in the integration kernels Eq. (5.11). Let us now take a closer look at the different elliptic topologies and the corresponding elliptic curves relevant for this publication. The topbox diagram has three elliptic sub-sectors corresponding to three different elliptic curves. They can be obtained from the maximal cuts in the Baikov representation [91][92][93]. We identify the three different elliptic curves by the labels a, b, and c according to the diagrams depicted in Fig. 9.
The elliptic curve E (a) is associated to the sunrise topology (Fig. 9a) and is probably the most well-known of the three. The corresponding z i in Eq. (5.19) are given as Iterated integrals involving only this topology can be cast in terms of iterated integrals of modular forms, and are well-suited for numerical evaluation [19,23,24,[94][95][96][97]. The elliptic curve E (b) is associated to the topbox sector itself (Fig. 9b) and the corresponding z i are given by The elliptic curve E (c) is associated to the remaining elliptic sub-sector which we dub the bubblebox sector (Fig. 9c). The corresponding quartic is defined by A more thorough analysis of the different sub-sectors and the corresponding elliptic curves can be found in [18].

Results for the new master integrals
Now that we have identified the various elliptic curves entering the computation, let us turn to the computation of the master integrals J 21 and J 22 , which are not known from the literature. As we have previously explained, the missing master integrals can be computed order by order in from the differential equation given in Eq. (5.9). We compute the master integrals iteratively in the expansion parameter by integrating the differential equation from a base point (x 0 , y 0 ) = (0, 1) (corresponding to s = ∞ and t = m 2 ). The integrals that need to be evaluated receive contributions from the elliptic topbox sub-sectors and hence contain elliptic iterated integrals themselves. The boundary values of the integrals at different orders in are inferred from the regular point x, y = (1, 1), where these two master integrals vanish. The integral J 21 is associated with elliptic curve a, and the first five orders in for J 21 read We have adopted the notation for the integration kernels appearing in the topbox family which was introduced in [18]. The integration kernels in Eq. (5.28) are given by We find that J 22 contains the elliptic sector 93 from [18] in one of its sub-topologies and is hence associated with the elliptic curve b. The first four orders in for J 21 can be written entirely in terms of MPLs and read The explicit dependence on the elliptic curve associated with the topbox sector enters for the first time at O( 4 ). The corresponding term involves 28 new integration kernels that have not appeared in the computations of the other master integrals and, accordingly, yields a rather large expression. We therefore refrain from showing it here explicitly and refer the interested reader to the ancillary files. Note that the differential 1-forms appearing in the iterated integrals are generally not path independent [85]. Nevertheless it is clear from a physical point of view that the master integrals do not depend on the particular choice of path chosen for the numerical evaluation of these iterated integrals. We have verified this path independence numerically using the methods explained in the following section.

Numerical evaluations and iterated integrals
In the previous sections, we discussed the analytic expressions of all the master integrals related to our process. In this section, we discuss how to numerically compute these integrals. The integrals containing only MPLs can be evaluated to a high precision, for example, using Ginac [98]. The numerical evaluation of iterated integrals associated with elliptic curves on the other hand is quite challenging. We explain how to compute these iterated integrals in two phase space regions, specifically, in the Euclidean and the physical region. Computing the iterated integrals in the Euclidean region in the small neighbourhood of the boundary point is relatively straightforward, as was also explained in [18]. For numerically evaluating these integrals in the physical region we need to analytically continue the associated elliptic curves across all the branch cuts appropriately. The analytic continuation of integrals having one parameter dependence has already been discussed in the literature [22,24,97]. However, analytically continuing all the integrals in our case, which includes the topbox integral associated to three different elliptic curves, is much more involved. With the objective of not digressing too much from the main goal of this paper, we discuss here only the ingredients essential for this computation.
The physical region for our case is governed by the following equations (m t = 1): where G is the Gram-determinant. In the coordinates x and y, these equations take the form  We integrate our iterated integrals from the base point (x = 0, y = 1) to any x and y over many segments and use path decomposition formula to evaluate the contribution over the whole path. The path segments need to be chosen while taking into consideration both spurious as well as physical types of singularities in the system.
The Euclidean region is defined as the region where the iterated integrals have only real contributions. The numerical evaluation of the iterated integrals that appear in our results for the master integrals is done for the Euclidean region as follows. We split the integration path from the boundary point (0, 1) to a generic kinematic point (x, y) in the Euclidean region into two pieces and use the path decomposition formula Eq. (5.14) to evaluate the integrals. It is better to use the coordinate system in Eq. 5.6 (x,ỹ) here.x rationalizes all the square roots associated with the non-elliptic kernels andỹ is used to bring the boundary point of integration to (0, 0). We choose the following paths: α, from (0, 0) to (x, 0) and β, from (x, 0) to (x,ỹ). For all the master integrals in our case, the integration along the first part gives only MPLs, which can be evaluated to high precision, as mentioned before. For the integration along the second part, we expand all integration kernels around y = 1, assumingỹ to be small. For integrating to a point y not close to 1, we may use more segments along y and use the path decomposition formula recursively.
To evaluate the master integrals in the physical region, we need to analytically continue the iterated integrals around the physical branch points. We use multiple (one-dimensional) path segments and use series expansion of the integrands on each of these paths to perform the integration. The choice of path is controlled by the radius of convergence of the series, which in turn depends on the singularities present in the kernels. It is important to choose properly the number of path segments and their sizes, so that the series solution converges properly on each of these segments and also the computation does not become too slow due to the presence of too many segments. Generally, the segments should not be larger than half the radius of convergence of the series expansion [20]. We again use the path decomposition formula to patch together the contribution from the different segments and compute the iterated integral over the whole path. For all the sunrise type kernels, there is a singularity while crossing the line y = 0 (ỹ = 1). For the elliptic kernels having both x and y dependence, there is a singularity on the line x = y, i.e., the line crossing over to the physical region. For the path shown in Fig.  11, both these type of singularities coincide. We therefore find it easier to integrate all the iterated integrals first along y and then along x. The kernels are no longer only MPLs on either of these paths but evaluate to iterated integrals over elliptic curves.
For all the integrals from the topbox family, we compute analytically in the physical region all the iterated integrals associated with the curves a and b, along with all the integrals depending on MPLs only. The periods of elliptic curve c are discontinuous on the line x = 0, unlike periods a and b. We reserve for a future publication a detailed explanation of the analytic continuation of the iterated integrals of the topbox family, which depend on multiple elliptic curves, together with the study of the analytic continuation of the integrals which depend also on the curve c.
Numerical integration to other regions in the phase space, not discussed above, can also be performed in a similar way. For example, integrating to the point (−0.5, 0.5) in Fig. 11 does not require us to cross any physical singularities if we integrate all the iterated integrals from the base point (0, 1). This can also be performed in the same way using the path decomposition formula recursively and does not involve a lot of complications. Another way to perform the integration in different phase space points is to start integrating the iterated integrals from a point already in that region. While this avoids the problem of analytic continuation it requires us to compute the boundary conditions at the new base points.

Functional Relations Among Iterated Integrals
In the previous sections, we have presented how we compute the planar master integrals contributing to the gg → tt scattering process in terms of MPLs and iterated integrals over elliptic curves. While this allows us to evaluate the ttgg helicity amplitudes, we find that there are potential redundancies in the functional basis of iterated integrals that might make the expression appear complex. Indeed, after integrating the differential equation, we found that the simple pole vanished numerically as expected but this cancellation was not reflected in the analytic expression. In order to achieve an explicit cancellation of all pole terms, some functional relations among the iterated integrals had to be applied. In this section we comment on functional relations among iterated integrals and the cancellation of the pole. For this purpose, we adopt again the notation for the integration kernels introduced in [18].
The functional relations we will discuss in this section are related to the appearance of derivatives of elliptic periods in the integration kernels. Take for example the kernel We see that the integration kernel a We can furthermore use integration by parts to derive other relations for iterated integrals involving the kernel a where f is an arbitrary integration kernel. Note that the first factor in the integral in Eq. (6.6) denotes a total differential and not the integration measure. If we can recast the product as a linear combination of integration kernels, we have successfully identified a relation among iterated integrals I which is closed under the minimal set of integration kernels we are using. We find the following equalities Plugging the above linear combinations into Eq. (6.7), we obtain three functional relations which we can use to simplify the pole terms. After doing so we are left with iterated integrals involving only rational kernels. These integrals evaluate to MPLs and cancel the ones already present in the expression, leading to the exact cancellation of the pole terms.
We can also use these relations to reduce the number of elements in the monomial basis of special function in the two-loop finite remainders with a closed top-quark loop A (2),N h . When mapping the master integrals appearing in A (2),N h to a basis of special functions, we found that there were 12025 monomials which, applying relations involving a (b) 3,3 as described above, was reduced to 11791 monomials. We notice that there are much fewer monomial basis elements appearing in the two-loop finite remainders F (2),N h . For example, for one of the subamplitudes with helicities + + +−, only 3586 monomial basis elements appear in the finite remainder with non-zero coefficients. Applying the relations discussed above we find that this number is reduced to 3158. Similar situations are also observed for the other subamplitudes and the + + ++ helicity configuration.
Note also that there are more kernels containing derivatives of elliptic periods. It is possible to write some of these kernels as a sum of a total derivative which captures all appearing derivatives of elliptic periods and a new kernel which does not contain the derivative of an elliptic period. By doing so, it might be possible to find further relations. Because of the large number of iterated integrals contributing to the scattering amplitude and the need to introduce new integration kernels, from our experience with a (b) 3,3 kernel, finding more relations between iterated integrals is necessary to simplify the finite remainder. We therefore reserve a thorough analysis of the possible functional relations for another publication.

Numerical Results
For benchmarking purposes, we present numerical results at a phase space point. We use a Euclidean phase space point that is considered in [18]  The rational representation of the momenta listed above is provided in the ancillary files.
The set of momenta given in Eq. We present results for the mass renormalised two-loop helicity amplitude for various fermion loop contributions specified in Eq. (2.4), normalised to the tree-level amplitudê with helicities λ i and f = 1, N l , N h , N 2 l , N l N h , N 2 h . Numerical results at the kinematic point given in Eq. (7.1) are displayed in Table 1. We use Ginac to evaluate MPLs that appear in the amplitude, while the numerical evaluation of iterated integrals involving elliptic kernels are done using the method described in Section 6.
We have performed several checks to validate the analytic results derived in this paper: • We compare numerical results in the Euclidean region for the squared matrix element obtained from the helicity amplitudes against the squared matrix element derived from an independent computation that directly computes the interference between the tree level and the two loop amplitudes in the HV scheme.
• We compare the finite remainder of the squared matrix element obtained from the helicity amplitudes against the results of [14] at the following physical point   with N l +N h active flavours in the running of α s has been used. In order to extract the finite remainder with N l +N h active flavours from [14], we subtract the IR poles in the CDR scheme from the renormalised amplitude given in Table 3 of [14]. As explained in Section 6, in this paper we have considered evaluation of the iterated integrals associated with the curves a and b in the physical region from analytical results while leaving the study of the analytic continuation of the integrals which depend also on the curve c for the future publication. Therefore, for the evaluation of the amplitude with a single closed top-loop contribution (A (2),N h ) in the physical region, we obtain numerical results of the mass renormalised amplitude in the master integral representation (Eq. (4.14)). In this representation, the master integrals containing elliptic curve c are evaluated using Fiesta [99] and pySecDec [100], while the remaining master integrals are evaluated from their analytical expressions. From this master integral representation we further remove the UV and IR poles numerically to obtain the finite remainder.

Conclusions
In this article we have computed a set of analytic helicity amplitudes for the planar corrections to top quark pair production via gluon fusion at two loops in QCD. While these amplitudes have been available for some time in alternative forms, this is the first time that amplitude level expressions using the massive spinor-helicity formalism have been obtained. We demonstrate that this form is suitable for processing with a rational parametrisation of the kinematics. This has the benefit that the algebra can be performed using finite field arithmetic and so combat the growth in algebraic complexity. This method will scale better than conventional approaches when considering additional final state particles.
Another important new ingredient is the inclusion of top quark loops using analytic expressions for the master integrals. Firstly we observed that the complete set of master integrals contributing to the amplitude were not available in the literature and computed them using the differential equation method. These integrals give rise to a set of special functions which have recently been presented in terms of iterated integrals over kernels including three elliptic curves. The growth in analytic complexity encountered in this sector is considerable and we presented a study of the numerical evaluation as well as analytic simplification leading to complete cancellation of the universal IR and UV poles.
We provide a complete set of independent finite remainders in the ancillary files accompanying the arXiv version of this article. We have also performed a cross check on the ingredients in our computation against the available numerical results.
Our work leaves a few open questions that motivate further investigation. The numerical evaluation of the iterated integrals is not at the same level of maturity as the commonly used MPLs. In particular the analytic continuation allowing a stable and efficient evaluation over the whole physical scattering region requires further study. We also notice that the iterated integrals obey additional functional relations. A complete understanding of these identities could and lead to a substantial reduction in the complexity of the final expressions and deserves additional study.

B Numerical results for one-loop amplitudes
We present in Table 2 numerical results for one-loop mass-renormalised amplitude at the kinematic point given in Eq.