One-loop QCD helicity amplitudes for pp → tt¯j\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t}j $$\end{document} to O(ε2)

We compute helicity amplitudes for the one-loop QCD corrections to top-quark pair production analytically in terms of a set of uniformly transcendental master integrals. We provide corrections up to O(ε2) in the dimensional regulator for the first time which are relevant at NNLO. Four independent pentagon integral topologies appear in the complete description of the colour structure for which we provide numerical solutions using canonical form differential equations and the method of generalised power series expansions. Analytic forms of the boundary values are obtained in all cases except one where we find a one-dimensional integral representation.


Introduction
Precision predictions for the production of a pair of top-quarks in association with a jet in hadron collisions is a high priority for current and future experimental measurements. Nextto-next-to-leading order (NNLO) corrections in Quantum-Chromodynamics (QCD) would allow percent level predictions for a wide variety of observables. The theoretical challenge and the degree of calculation complexity for such predictions remains extremely high.

JHEP06(2022)066
Predictions for top quark pair production in association with multiple jets [12] or matched to a parton shower [13][14][15] have been made possible thanks to the latest generation of automated numerical tools. The pp → ttj process is of particular interest since it is extremely sensitive to the top quark mass [16,17].
Precision predictions at NNLO are currently only available for the four particle process pp → tt. Advanced techniques for the subtraction of infrared divergences [18] have enabled a comprehensive range of phenomenological studies [19][20][21]. The amplitude level ingredients for these predictions are largely known analytically [22][23][24][25][26][27][28][29][30] although there are still a small number of non-planar double-virtual contributions that are only known numerically.
In this article we present one previously missing ingredient relevant for a next-to-nextto-leading computation of pp → ttj: the expansion of the one-loop helicity amplitudes up to O(ε 2 ) in the dimensional regulator. This requires the computation of new pentagon integrals that first appear at O(ε) and are one of the new results presented here.
Helicity amplitudes (including decay information for the top-quark pair in the narrow width approximation) for this process have not been presented analytically before. Oneloop expressions for pp → tt production were computed in this formalism using unitarity based methods and led to relatively compact expressions [31]. The motivation to do so here for the high multiplicity process is to get a sense of the complexity that might arise in an analytic two-loop computation of pp → ttj. The new loop integrals appearing at O(ε) depend on genuine five-point kinematics for the first time. While at one-loop all the special functions are of a polylogarithmic form, the alphabet is quite complex and efficient evaluation and analytic continuation to physical kinematics is challenging. In this context recent progress has been made to compute analytically five-point one-loop integrals with massive external legs and internal propagators [32]. In this article we explore the technique of generalised series expansions [33], as implemented in the software DiffExp [34], to numerically solve the differential equations for the master integrals. Such a technique would be applicable to two-loop integrals even in the presence of non-polylogarithmic forms, and it has been exploited recently for several processes [35][36][37][38]. This approach as well as related methods for the numerical solution of differential equations for master integrals [39][40][41][42][43] have been of particular interest recently due to their wide range of applicability.
The rational coefficients of the special functions also represent a step up in analytic complexity in comparison with previously considered two-loop five-point amplitudes . We present a complete set of partial colour amplitudes in terms of master integrals valid to all orders in the dimensional regulator. These objects are considerably more complex than the four-dimensional limits and we employ a finite field reconstruction technique [66,67] and a rational parametrisation of the kinematics based on momentum twistors [68] to overcome the algebraic complexity.
Our paper is organised as follows. We begin by reviewing the colour decomposition of the amplitudes in both ttggg and ttqqg channels and describe the infrared and ultraviolet pole structure. We then describe the finite field reconstruction approach taken to extract the independent helicity amplitudes. We then turn our attention to the evaluation of the master integrals. We present canonical form differential equations for the four independent topologies appearing in our process. The computation of the boundary terms is described JHEP06(2022)066 and the numerical evaluation using generalised series expansions with DiffExp is presented. Finally we present some numerical results before giving an outlook for the future.

Colour decomposition and infrared pole structure
We choose to define the two partonic channels for pp → ttj with all momenta out-going. Evaluation for physical kinematics can be performed using the appropriate analytic continuation. We write the amplitudes according to the colour decomposition [69]. Therefore, for the process 0 →ttggg we have: Here we have used g s to denote strong coupling and taken an overall normalisation are the partial amplitudes which appear in the full amplitude as sums over permutations of the momenta. S 3 indicates the six permutations of the three gluons while S 3 /Z 2 and S 3 /Z 3 are smaller symmetry groups with 3 and 2 elements respectively. The SU(N c ) colour structures are written using the fundamental generators (t a )j i where a = 1, · · · , 8 are indices of the adjoint representation, while i,j = 1, 2, 3 are indices in the fundamental and anti-fundamental representation respectively.
Following the same conventions, we colour decompose the process 0 →ttqqg as (see for example [70]), Each of the partial amplitudes is further decomposed into a polynomial in N c and the number of light and heavy flavours, n f and n h = 1 respectively. Suppressing the momentum JHEP06(2022)066 arguments we have for the 0 →ttggg channel, and for 0 →ttqqg channel where X = 1, · · · , 4. The kinematics for these processes is: where p 1 and p 2 are the momenta of the external top quarks, p 3 and p 4 are the momenta associated either to the a pair of gluons or a pair of massless quarks, and p 5 is the momentum of the remaining gluon. All the particles are on-shell and the top quarks are considered to be massive, with m t the top mass. Finally, throughout this paper, we work both with the kinematic invariants d ij and s ij as defined in (2.12).

Infrared singularities
Catani, Dittmaier and Trocsanyi (CDT) were the first to present a closed formula for the universal infrared (and ultraviolet) pole structure of an arbitrary one-loop amplitude with massless and massive QCD partons [71]. Using the colour space notation [72] the factorisation of the infrared poles can be denoted simply as, For renormalised amplitudes the pole operator I n is defined as, 1 (2.14) 1 We omit the imaginary parts in our reproduction of the In operator. For a correct treatment across the full physical phase-space the prescription is given in ref. [71]. At the test points we provide, the form given and Mathematica's internal prescription are sufficient to find agreement up to O(ε −1 ).

JHEP06(2022)066
where we have followed the normalisation conventions from eq. (2.2). 2 The function V ij arises from the soft singularities which contains colour correlations of the form T i · T j , i and j are massive comes from the kinematic threshold for the production of a top quark pair. Since we only have two massive partons in our process with the same mass we have a single threshold β 12 The finite parts of the function V ij will not play a role in the cross checks of our computation though they are important to ensure the correct small mass limits, as stated in [71]. The functions Γ j arise from the hard collinear region and depend on the anomalous dimensions of the partons. Our amplitudes are computed including wave-function renormalisation but excluding coupling renormalisation and therefore additional UV poles proportional to the QCD β function are present in the expressions.
We will compute all partial amplitudes in terms of master integrals valid to all orders in ε. The verification of the infrared pole structure is therefore an extremely strong check on the validity of our expressions. The inclusion of the wave-function renormalisation counterterms ensures the amplitude is gauge invariant which also provides a strong cross check.
Explicit evaluations of the CDT formula, including for the two partonic channels, into the partial decompositions eq. (2.1) and eq. (2.3) are given in appendix B.

Helicity amplitude setup
In this section we describe the computational set up for the helicity amplitudes. This makes use of the well known spinor-helicity formalism for massless and massive particles.

Helicity amplitudes
The helicity states for the massive fermions are computed using the standard decomposition along an arbitrary reference direction [73]: where p is the massive fermion momentum, m is the mass, n is the arbitrary reference direction and p = p − m 2 2p·n n. Since the direction n is arbitrary the positive helicity state

JHEP06(2022)066
is related to the negative helicity state through the transformation n ↔ p together with a normalisation factor accounting for the change in spinor phase. For further details of this relation, and other aspects of the massive spinor-helicity formalism used in this article, we point the reader to [30] and references therein.
We perform an analytic reconstruction in the minimal set of six on-shell variables making use of the following basis for spin structures, The decomposition involves a phase factor Φ to account for the massless parton helicities, four basis functions Θ i for the spin dependence of the top-quark pair and the associated subamplitudes A . This representation and notation has been introduced in the recent study of top-quark production [30]. The functions Θ contain all dependence on the arbitrary reference vectors introduced to define the positive helicity massive fermions. Such a spin basis is not unique and the normalisations for the Θ functions have been chosen such that all subamplitudes have the same dimension and are free of any spinor phase. This form is sufficient to account for the decays of the top quarks in the narrow width approximation [3,74].
For the amplitudes considered in this article the explicit forms for Φ and Θ are: and, The amplitudes with a massless fermion pair use the same spin decomposition (3.9) but

Rational phase space parametrisation
Our computation uses a rational phase-space parametrisation together with a numerical sampling of the relevant set of Feynman diagrams using modular arithmetic. The generation of this rational parametrisation uses the momentum twistor formalism [68], a technique that has been applied numerous times in similar amplitude computations. For the case of a top-quark pair plus three massless partons the method is essentially the same as the one described in ref. [30] for a top-quark pair plus two massless partons. We begin by generating a rational parametrisation for configurations of seven massless particles (11 free variables) with momenta q 1 , · · · q 7 . 3 The five particle system for tt plus three partons can then be written, with additional constraints to ensure the massive momenta p 1 and p 2 are on-shell. Specifically we solve the constraints: Having found the rational parametrisation we change variables to, (3.14) In the last variable we have introduced the notation p ij = p i + p j . We note that the only dimensionful variable s 34 can be set to 1 and restored easily through dimensional analysis. It is not possible to use the top quark mass as a variable without introducing square roots and hence we choose a spinorial trace. For completeness we present explicitly the map

Reduction to master integrals
Our amplitude computation strategy follows the method applied to recent two-loop computations [30]. The colour-ordered helicity amplitudes are first generated from Feynman diagrams using QGraf [78] which are then processed using a combination of Mathematica and FORM [79,80] scripts. The Spinney package is used to process parts of the numerator algebra [81]. The integral topologies are identified from the loop momentum dependent propagators and the coefficients of the loop-dependent numerators are computed using the momentum twistor parametrisation. A tensor integral representation for each diagram is obtained using a basis of irreducible numerators for the maximal cut topologies.
Since the maximal cut has four independent momenta, no transverse integration [82] step is required in contrast to the method of reference [30]. The diagram numerators are computed using a symbolic value for the spin dimension d s = g µ µ and our results are presented as an expansion around d s = 2, (3.26) The dimension d s used in the numerator algebra is kept separate from the loop integration dimension d = 4 − 2ε in order to account for the scheme dependence between the FDH (d s = 4) and tHV (d s = 4 − 2ε) schemes. The diagram numerators are then reduced to master integrals via integration-by-parts identities [83,84] following the Laporta algorithm [85]. The complete reduction is implemented in LiteRed [86] and FiniteFlow [67] so that numerical sampling with modular arithmetic can be used to reconstruct the coefficients of the master integrals directly without analytic intermediate steps. Wave-function renormalisation must also be performed in order to obtain a gauge invariant result. We generate diagrams with the counter-term insertions which are added to the one-loop numerators as shown in figure 1. The counter-terms are written in terms of loop integrals JHEP06(2022)066 valid to all orders in ε. The Feynman rule for the counter-term insertion can be written as, The right hand side of this equation consists of: a line representing the Feynman rule for a massive fermion propagator, a rational function of N c , d s and ε, and a wave-function bubble graph representing the Feynman integral with one massive and one massless propagator and a mass scale of m 2 t . We note that using this integral form for the counter-term, the amplitude is already guage invariant at the level of master integrals. This ensures that only gauge invariant quantities are reconstructed analytically and allows us to sidestep the issue of including external leg corrections with on-shell ingredients [8,87,88].
The coefficients of the master integrals are functions of the dimensional regularisation parameter ε and the six free parameters in the rational phase-space. Since this is a oneloop problem the evaluation time for the amplitude within the FiniteFlow setup is quite fast and so functional reconstruction of high degree polynomials is feasible. Nevertheless, we find that rationalising the phase-space increases the polynomial degree significantly so in addition we apply linear relations and a univariate partial fractioning to the master integral coefficients before reconstruction as has been effective in many cases with massless propagators(for example [50]). We apply the univariate reconstruction method and the algorithm for linear relations described in ref. [62]. The first step in this method requires a matching of denominator (and numerator) factors for which we build an ansatz from a set of spinor products, Gram determinants and other denominators appearing in the differential equations which we describe in the next section. We have used the following JHEP06(2022)066 kinematic structures to generate our factor ansatz, where we have introduced a notation for the 3-mass triangle Gram determinants, ∆ 3 and the Cayley matrix associated with the pentagon integral with four internal masses, We note that in spinor-helicity variables, many of the Gram and Cayley determinants factorise and so we don't need to specify all Gram and Cayley determinants explicitly.
To generate an ansatz that matches all denominators this list is permuted over an (overcomplete) set of six permutations of 3, 4, 5 and two permutations of 1, 2 (12 in total) and duplicate entries are removed. To match the polynomial factors this list is evaluated using the rational momenta parametrisations from which a list of independent polynomials can be extracted.
We summarise our reduction strategy as follows: we apply the rational kinematic parametrisation to the processed Feynman graphs for the four different spin projections. Each of these projected amplitudes are reduced to master integrals and reconstructed using FiniteFlow. After reconstruction, the projected amplitudes are used to construct the sub-amplitudes, again with the help of reconstruction over finite fields, linear relations and univariate partial fractioning.

Computation of the master integrals
There are four distinct pentagon function topologies appearing in the amplitudes as shown in figure 2. To find the minimal set of master integrals (MIs) which describes each topology we perform Integration-By-Parts (IBPs) reduction [84,89], as implemented in the software LiteRed [86,90] and FiniteFlow [67]. We find that the topologies T 1 , T 2 , T 3 and T 4 are described, respectively, by 15, 21, 17 and 19 MIs (see figures 3, 4, 5 and 6). By employing symmetry relations among the different topologies, and their permutations, we find that JHEP06(2022)066 the amplitudes can be written in terms of minimal set of 130 MIs across the four topologies. The evaluation of the MIs that appear in the amplitudes is discussed in section 5.1.
We compute the MIs, f (x, ε), by means of the differential equations method [91][92][93]. Specifically, we work with a system of differential equations in canonical form [94]: where d is the total differential with respect to the kinematic invariants, The matrix A( x) is a linear combination of logarithms: where c i are matrices of rational numbers and α i ( x) are algebraic functions of the kinematic invariants x.
A major feature of the systems of differential equations for the four pentagon topologies is that they depend on the following set of square roots: where the argument a 1 can be functions of the kinematic invariants, P and Q are momenta and G ij ( v) = 2v i · v j is the Gram matrix.

JHEP06(2022)066
We choose to solve the systems of differential equations for the MIs using the generalized power series expansion method [33], as implemented in the software DiffExp [34]. Although this method furnishes a semi-analytic solution to the MIs, it has the advantage of allowing a fast and high precision numerical evaluation. Moreover, the analytic continuation of the solution is easier with respect to an analytic approach. Indeed, while it could be possible to linearize the square roots system 4 (4.4) with a transformation of the kinematic invariants, and obtain an analytic solution in terms of polylogarithmic functions (MPLs) [95,96], the system of differential equations will involve polynomials of high degree in the linearized variables. This feature impacts significantly the computation since the system of differential equations in the new set of variables is too large to be handled efficiently. In addition, the determination of the phase-space regions, and therefore the analytic continuation, is more complicated.
We finish the first part of this section by discussing a few more details of our computation. Firstly, we would like to stress that we reconstruct the systems of differential equations for the four pentagon topologies exploiting finite fields methods, implemented in FiniteFlow, for a basis of master integrals, f that does not contain the square roots (4.4): (4.5) In doing so, we obtain a system of differential equations that is not in a canonical form, but we also avoid dealing with square roots in the reconstruction procedure. Then, an ε-factorized form can be achieved by a rotation of the basis of MIs, f = B( x) f , under which the matrix d A ( x, ε) transforms as: where B( x) is a diagonal matrix whose entries are the square roots (4.4). Obtaining the matrix A ( x, ε) can in general be complicated for multi-loop cases but for this one-loop case it is straightforward. The rotation matrix B( x) is easily obtained using information from the maximal cuts of the topologies. Finally, we comment on the computation of the boundary conditions for the system (4.1). We compute the boundary conditions, for a minimal subset of the MIs, by direct integration of their Feynman parameter representation at the kinematic point: This step is performed using the linear-reducibility strategy as implemented in Hyper-Int [97] with the help of PolyLogTools [98]. High precision numerical boundary values are obtained for the remaining integrals using DiffExp. We detail the boundary condition computation for each topology in the following subsections. We include supplementary material containing the systems of differential equations

The T 1 topology with one massive internal propagator
There are 15 master integrals in the topology T 1 as shown in figure 3. The integrals are defined as: where (4.10) The basis of canonical MIs is chosen to be: 1,1,1,1,1 , 1,0,1,1,1 , 1,1,0,1,1 , 1,1,0,1,0 , We compute analytically the boundary conditions for all the MIs in this topology, for which we obtain expressions in terms of MPL functions. Then we use GiNaC to evaluate them numerically with high precision (100 digits). The integral f T 1 15 is equivalent to a tadpole, and therefore we can write its boundary value exactly:

The T 2 topology with 4 massive internal propagators
Topology T 2 is described by 21 master integrals as shown in figure 4. The integrals are defined as:

JHEP06(2022)066
where (4.14) The canonical basis for the topology T 2 is chosen to be: As for topology T 1 , we compute analytically the boundary conditions for all the MIs in topology T 2 but for the pentagon f T 2 1 , for which we obtain an expression in terms of JHEP06(2022)066 one-parameter integrals. Moreover, since we do not have to perform any new computation for these integrals.

The T 3 topology with 2 massive internal propagators
Topology T 3 is described by 17 master integrals as shown in figure 5. The integrals are defined as: a 2 ,a 3 ,a 4 ,

JHEP06(2022)066
where The canonical basis for the topology T 3 is chosen to be: 1,0,1,1,1 , 1,1,0,1,1 , 2,0,0,1,0 , 0,2,0,0,1 , The boundary conditions for the bubble integrals f T 3 17 , f T 3 16 , f T 3 15 , f T 3 13 are known from topology T 1 , while the values for the integrals f T 3 14 , f T 3 12 can be obtained numerically using Dif-fExp. For example, the boundary condition for f T 3 14 can be obtained from f T 3 13 by evolving, using DiffExp, its value at x 0 to a point x σ 0 , whose value is determined by an appropriate permutation of the external momenta. Indeed, f T 3 14 is the same integral as f T 3 13 just in a different channel, as it can be seen from figure 5. We also exploit this strategy for other integrals in T 3 and T 4 .

JHEP06(2022)066
In T 3 , only the boundary values for the integrals f T 3 1 , f T 3 2 , f T 3 6 and f T 3 9 need to be computed explicitly through direct integration. All other boundary values have already been computed in previous topologies or through numerical evaluation using DiffExp. In particular, the integrals f T 3 8 and f T 3 4 appear already in other topologies, 20) and the boundary values for the remaining MIs are evaluated with DiffExp in the following way: • the boundary value for f T 3 11 is obtained from f T 2 15 ; • the boundary value for f T 3 10 is obtained from f T 2 9 ; • the boundary value for f T 3 7 is obtained from f T 1 8 ; • the boundary value for f T 3 5 is obtained from f T 3 2 ; • the boundary value for f T 3 3 is obtained from f T 1 4 .

The T 4 topology with 3 massive internal propagators
Topology T 4 is described by 19 master integrals as shown in figure 6. The integrals are defined as: where The canonical basis for the topology T 4 is chosen to be: 2,1,0,0,0 .

Numerical results for the master integrals
In this section we discuss our results for the numerical evaluation of the MIs performed with DiffExp. The amplitudes depend on 130 independent MIs across all the four pentagon topologies, and their permutations. Instead of evaluating the whole system of 130 MIs at once we evaluate each topology separately. Since the number of MIs inside each topology is at most 21, this approach allows us a faster numerical evaluation, as we can evaluate in parallel all the topologies and their permutations. The timing to get numerical values for all the topologies and permutations is within a range of ∼ 30 minutes to ∼ 1 hour for phase-space point, on a laptop, requiring an accuracy of 16 digits. We stress that for phenomenological applications the performances can be improved by building a precomputed grid of points as boundary values [35,36]. Given a point x a , we can evaluate all the MIs, and hence the amplitude, at that point as follows. The standard ordering of the external momenta for topologies T 1 and T 2 is (1, 2, 3, 4, 5), while for topologies T 3 and T 4 it is (1, 3, 2, 4, 5). The permutations of these topologies are given by all the possible permutations of the momenta (3,4,5). Therefore, evaluating the MIs that belong to the permutations of T 1 ,T 2 ,T 3 and T 4 is equivalent to evaluating the MIs that describe the topologies in the standard orderings (1, 2, 3, 4, 5) and (1, 3, 2, 4, 5) at a kinematic point, x σ a , which is given by the corresponding permutations of the kinematic invariants.
In order to clarify this procedure we discuss the following example. We consider the permutation (1, 2, 4, 3, 5) for the topology T 1 . The permutation of the external momenta:

JHEP06(2022)066
implies the following transformation for the kinematics invariants: This means that evaluating the permutation (1, 2, 4, 3, 5) of T 1 at the point: is equivalent to evaluating T 1 , in the standard ordering (1,2,3,4,5), at the point: This procedure allows us to evaluate all the permutations of a given topology starting from the system of differential equations for the topology in the standard ordering. Moreover, this strategy has also been used to compute the boundary conditions for some MIs as discussed in the previous section.
In order to verify the correctness of our computation we performed different checks comparing our results against numerical values for the MIs obtained by means of sector decomposition techniques, as implemented in the software PySecDec [99].

Amplitude results
The explicit analytic forms for the partial helicity amplitudes, broken into subamplitudes according to eq. (3.2), are provided in the supplementary material. Due to the large overall size we do not attempt to provide any typeset expressions in the paper. The coefficients appearing in the subamplitudes have been collected and linear relations between them determined. To provide a relatively compact format common factors in the linearly indpendent rational coefficients are identified and presented as a set of replacement rules.
The univariate partial fractioning in the variable x 5123 was quite effective in reducing the total number of sample points required in the reconstruction. The maximum total degree appearing in the most complicated sub-leading colour amplitudes was O(100) before partial fractioning and linear relations. This reduced to O (20) in the final expressions. Nevertheless, the appearance of high degree polynomials does indicate that extensions to higher loops will be challenging since the fast evaluation of the one-loop input enables to handle such expressions without restrictions on computational resources.
The supplementary material also provide two example scripts demonstrating the evaluation of the amplitudes using the numerical results for the master integrals and validation of the universal pole structure. All evaluations have been performed in Mathematica where we can obtain O(100) accurate digits without issues. Since the main use case of the new O(ε) and O(ε 2 ) terms will be in the subtraction of divergences in two-loop amplitudes, we have not attempted to provide an efficient evaluation of the amplitudes for use at NLO.

JHEP06(2022)066 6 Conclusions
In this article we have presented a computation of all one-loop helicity amplitudes of the process pp → ttj evaluated to O(ε 2 ). The expansion to higher order in ε allows us to look at the complexity of the NNLO terms for the first time. Applying finite field reconstruction techniques demonstrates that the algebraic complexity of this problem may be within reach. The analytic complexity coming from the loop integrals was easily overcome using the combination of canonical form differential equations and subsequent evaluation using generalised series expansions in DiffExp. The boundary terms are provided in analytic form up to weight four using MPLs except for the pentagon master integral with four internal masses which is presented as a one-paramater integral.
It will be interesting to see how automated approaches to loop integral evaluation using the numerical evaluation of the differential equations develop. Since mathematical bottlenecks in the understanding of elliptic structures and the difficulties of dealing with long and complicated alphabets may be sidestepped, the method has substantial advantages over fully analytic approaches. Nevertheless this comes at the cost of numerical performance and the determination of the boundary values will still be a major issue. Recent attempts to automate the evaluation of boundary terms using sector decomposition have been successful [100] although the numerical accuracy is probably not yet sufficient for the full phase space.
There are clearly important issues that should be addressed in order to overcome challenges at two-loops. Amplitudes in d = 4 − 2ε dimensions are substantially more complicated than their four-dimensional limits. The identification of an analytic function basis such that expansion in ε can be taken and the subtraction of poles can be performed is likely to be an essential ingredient. We also observe a high degree of algebraic complexity stemming from the global choice of rational kinematic parametrisation. This is particularly evident in the sub-leading colour partial amplitudes in which many permutations of the master integral topologies appear.
Despite significant challenges ahead, the work presented here motivates further investigation into analytic or semi-analytic approaches to high precision pp → ttj amplitudes and cross-sections.

A Generalized series expansion method: a brief review
For the reader's convenience, we start this section with a short review of the method of generalized power series [33,35], implemented in [34], which has been exploited to evaluate the MIs numerically.
The generalized series expansion method allows us to evaluate the solution to the system (4.1), at a point x a , from the knowledge of the solution at some boundary point x 0 . This is done is the following three steps: • Step 1: we split the integration path into segments; • Step 2: we find a solution inside each segment by expanding in series the system of differential equations; • Step 3: we evaluate the solution in the point x a connecting the local solutions of the path γ(t).
The solution to the canonical system (4.1) can be written as a series expansion in ε: where: and we assume that the solution is described by some variable t which parametrizes the path, γ(t), which connects the points x 0 and x a : As already mentioned, the first step consists in splitting the path γ(t) into segments where {t i } is the set of points in which we are going to expand the system of differential equations, and r i is the radius of convergence of the series inside each segment. The segments S i can be chosen from the knowledge of the singular points of the differential equations. In particular we can have both real: and complex-valued singular points: Therefore, we can choose the expansion points to belong to the set R ∪ C r , where C r is a set of regular points:

JHEP06(2022)066
and the radius of convergence, r i , can be defined as the distance of t i to the closest element t p , with p = i.
In the second step of the method we determine local solutions to the differential equations inside each segment S i . This is done by expanding the system of differential equations around the point t i : where A l are constant matrices. Then, exploiting the general expression (A.2) for the k-weight of the solution, we obtain: for the local solution, f (k) i (t), inside the segment S i . We point out that working with a system of differential equations in canonical form implies that the integrals that appear in (A.8) are of the form: Consequently, (A.8) is given by the expression: where the matrices c (i,l 1 ,l 2 ) k depend on the boundary conditions for the system and on the constant matrices A l in (A.7).
Finally, as last step of the procedure, the global solution on the path γ(t) can then be approximated as: where N is the total number of segments, and f (k) i (t) is the k-weight of the local solution, inside the segment S i , written as a truncated series expansion, around some point t i , with radius of convergence r i .

A.1 Analytic continuation
In this part we briefly discuss the analytic continuation within the framework of the generalized series expansion method and, in particular, in its DiffExp implementation.
Analytic continuation has to be performed when a singularity of the differential equations matrix, dA( x), is crossed along the integration from the boundary point, x 0 , to the evaluation point x a . We encounter two kinds of singularities:

JHEP06(2022)066
• Case I: Logarithmic singularities, this type of singularities arise from simple poles in the system of differential equations; • Case II: Square roots singularities, this type of singularities appear as the canonical basis of MIs involves square roots of the kinematic invariants.
In order to determine all possible logarithmic singularities, we consider the system of differential equations for a basis of MIs without square roots normalization involved, then we perform a multivariate partial fraction on the system exploiting the Mathematica package MultivariateApart [101]. In doing so, we obtain a set of irreducible polynomials in the kinematic invariants, P( x), which describes the simple poles structure of the system (A.7). Instead, the square roots singularities, S( x), are given by set of square roots (4.4) that define the canonical basis of MIs. Therefore, the full set of singularities for the analytic continuation is given by the set of polynomials P( x) ∪ S( x).
Once the full set of singularities is known we perform the analytic continuation as follows. We assign a small imaginary part to the kinematic invariants and the top mass: then we substitute (A.12) into the polynomials P( x) ∪ S( x), we expand with respect to δ and we keep just the linear term. As an example of this procedure we consider the logarithmic singularity As a consequence, once we substitute (A.14) into (A.13), we obtain that the singularity (A.13) is analytic continued as:

B Explicit form of the infrared poles for the partial amplitudes
In this appendix we give the explict form for the Catani-Dittmaier-Trocsanyi formula eq. (2.14) for the partial colour amplitudes. We use the following short hand for the logarithms that appear,