Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel
Simon Badger, Ekta Chaubey, Heribertus Bayu Hartanto, Robin Marzucca
PPrepared for submission to JHEP
Two-loop leading colour QCD helicity amplitudesfor top quark pair production in the gluon fusionchannel
Simon Badger, a Ekta Chaubey, a Heribertus Bayu Hartanto, b Robin Marzucca c a Dipartimento di Fisica and Arnold-Regge Center, Università di Torino, and INFN, Sezione diTorino, Via P. Giuria 1, I-10125 Torino, Italy b Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom c Institute for Particle Physics Phenomenology, Department of Physics, Durham University, DurhamDH1 3LE, United Kingdom
E-mail: [email protected], [email protected],[email protected], [email protected]
Abstract:
We present a complete set of analytic helicity amplitudes for top quark pairproduction via gluon fusion at two-loops in QCD. For the first time, we include correctionsdue to massive fermion loops which give rise to integrals over elliptic curves. We presentthe results of the missing master integrals needed to compute the amplitude and obtainan analytic form for the finite remainders in terms of iterated integrals using rationalisedkinematics and finite field sampling. We also study the numerical evaluation of the iteratedintegrals. a r X i v : . [ h e p - ph ] F e b ontents ¯ ttgg amplitudes 33 Helicity amplitudes for massive fermions 5 Precise understanding of the top quark in the Standard Model is one of the highest prioritiesin modern collider experiments. The large mass of the top quark makes it a dominantingredient in understanding the fundamental forces and in particular of the electroweaksymmetry breaking mechanism. A precise determination of the top quark mass is imperativedue to the sensitivity of the Higgs potential on its value.The QCD corrections to top-quark pair production have a long history going backto the first total cross section computation at NLO for on-shell top-quarks computed byNason, Dawson and Ellis [1, 2]. Thanks to a monumental effort on the part of the theo-retical community, fully differential NNLO predictions for top anti-top pair production athadron colliders are now available for comparisons with the latest experimental data [3–7].These intensive computations require two-loop → , one-loop → and tree-level → – 1 –artonic amplitudes to be combined using a consistent UV and IR subtraction scheme toobtain a finite result. The first complete computation used the ‘STRIPPER’ sector decom-position method to subtract IR divergences and has enabled the most complete theoreticaldescriptions of top-quark dynamics [8]. Very recently a second independent computationhas been obtained using the q T subtraction method [9–12].Up to now the only complete two-loop amplitudes have been available numerically [13–15]. The main reason for this is that a more complicated class of special functions begin toappear in amplitudes which contain internal masses. These functions have been identifiedto involve integrals over elliptic curves [16–24], and lie on the boundary of our currentmathematical understanding. Over the last few years, motivated by the increasing demandof precise theoretical descriptions including mass dependence, there has been substantialprogress 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 dependon elliptic sectors and there has been substantial effort to find a complete analytic form forthe squared amplitudes [25–32].In this paper we present a set of helicity amplitudes for top-quark pair production in theleading colour approximation. For the first time we also include contributions from heavyfermion 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 diagramswith finite field arithmetic [33–38], we also study the numerical evaluation of the finalamplitudes. The helicity amplitudes presented contain complete information about topquark decays in the narrow width approximation. This helicity amplitude technique hasbeen used successfully at one-loop [39–42]. Spin correlations can also be treated efficientlyusing the spin density matrix approach [15, 43, 44]As is always the case, computations at this perturbative order contain a large number ofsteps, each usually with some technical bottlenecks to overcome. The analytic computationsof the planar master integrals, including those with an internal massive loop, are availablein the literature for almost every topology. Yet, for the construction of the amplitude it wasnecessary to include one final two-loop integral topology which contained two elliptic masterintegrals. The (canonical form) differential equation [21, 45–49] and solution in terms ofiterated integrals was obtained and combined with the integration-by-parts reduced helicityamplitudes.Throughout the computation we made use of finite field arithmetic to find an efficientsolution to the system of integration-by-parts identities. These techniques have shown tobe particularly efficient for massless amplitudes with many external scales. In this paperwe 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 twistorsand 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 ofon-shell, gauge invariant sub-amplitudes that can be computed using on-shell top-quarkkinematics. We then check these results against previous computations for the interferencewith tree-level diagrams in the conventional dimensional regularisation (CDR) scheme.– 2 – (1) , A (1) ,N l A (1) ,N h Figure 1 : Sample Feynman diagrams corresponding to various internal flavour contribu-tions at one loop as specified in Eq. (2.3). Red lines, black spiral lines and black linesrepresent massless quarks, gluons and top quarks, respectively. A (2) , A (2) ,N l A (2) ,N h A (2) ,N l A (2) ,N l N h A (2) ,N h Figure 2 : Sample Feynman diagrams corresponding to various internal flavour contribu-tions at two loops as specified in Eq. (2.4). Red lines, black spiral lines and black linesrepresent massless quarks, gluons and top quarks, respectively. ¯ ttgg amplitudes We consider a scattering process involving a pair of top quarks and two gluons → ¯ t ( p ) + t ( p ) + g ( p ) + g ( p ) , where p = p = m t and p = p = 0 . The kinematic invariants for this process are thetop-quark mass m t , and the two Mandelstam variables s = ( p + p ) , t = ( p + p ) . (2.1)In this work we consider the leading colour contributions of the ¯ ttgg amplitude up totwo-loop level, where at two loops, only planar configurations arise. The colour decompo-sition of the leading colour L -loop ¯ ttgg amplitude is given by A ( L ) (1 ¯ t , t , g , g ) = n L g s (cid:20) ( T a T a ) ¯ i i A ( L ) (1 ¯ t , t , g , g ) + (3 ↔ (cid:21) , (2.2)where n = m (cid:15) α s / (4 π ) , α s = g s / (4 π ) , m (cid:15) = i (cid:0) π/m t (cid:1) (cid:15) e − (cid:15)γ E , g s is the strong couplingconstant and ( T a ) ¯ ji are the fundamental generators of SU ( N c ) .– 3 – (0)mct A (1) , A (1) ,N h mct A (0)2mct Figure 3 : Sample Feynman diagrams with mass counterterm insertions at tree level andone loop that contribute in Eqs. (2.5) and (2.6). Black spiral lines and black lines rep-resent gluons and top quarks, respectively. Circled crosses indicate the mass counterterminsertions.The partial amplitudes can further be decomposed according to their internal flavourstructure, A (1) (1 ¯ t , t , g , g ) = N c A (1) , + N l A (1) ,N l + N h A (1) ,N h , (2.3) A (2) (1 ¯ t , t , g , g ) = N c A (2) , + N c N l A (2) ,N l + N c N h A (2) ,N h + N l A (2) ,N l + N l N h A (2) ,N l N h + N h A (2) ,N h . (2.4)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 inFigs. 1 and 2.We first define the gauge invariant L -loop mass-renormalised amplitude A ( L )mren , wherethe top-quark mass renormalisation is included A (1)mren = A (1) − δZ (1) m A (0)mct , (2.5) A (2)mren = A (2) − δZ (2) m A (0)mct − δZ (1) m A (1)mct + (cid:0) δZ (1) m (cid:1) A (0)2mct . (2.6) δZ ( L ) m is the L -loop mass renormalisation counterterm, A (0)mct ( A (1)mct ) represents an ampli-tude 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 masscounterterm insertions to top-quark propagators. The mass counterterm insertion proce-dure leads to the presence of double and triple top-quark propagators in the tree levelamplitude and double top-quark propagators in the one-loop amplitude. Sample diagramswith mass counterterm insertions at tree level and one loop are displayed in Fig. 3.The fully renormalised ¯ ttgg amplitude can be obtained by including renormalisationsof the top-quark and gluon wavefunctions as well as the strong coupling α s A (1)ren = A (1)mren + (cid:0) δZ (1) t + δZ (1) g + δZ (1) α s (cid:1) A (0) , (2.7) A (2)ren = A (2)mren + (cid:0) δZ (1) t + δZ (1) g + 2 δZ (1) α s (cid:1) A (1)mren + (cid:0) δZ (2) t + δZ (2) g + δZ (2) α s + δZ (1) t δZ (1) g + 2 δZ (1) t δZ (1) α s + 2 δZ (1) g δZ (1) α s (cid:1) A (0) . (2.8)– 4 – Z ( L ) t , δZ ( L ) g and δZ ( L ) α s are the L -loop top-quark wavefunction, gluon wavefunction andstrong coupling renormalisation constants, respectively. All the renormalisation constantsrelevant for our calculation are given in Appendix A.The renormalisation procedure removes the ultraviolet (UV) singularities present inloop amplitudes. The remaining divergences that are of infrared (IR) nature can be pre-dicted from universal IR behaviour of QCD amplitudes [50–52], which, in the presence oftop quarks, needs to be extended to massive case [53–55]. The UV renormalised amplitudeat one and two loops can be divided into singular and finite terms |A (1) n (cid:105) = Z (1) |A (0) n (cid:105) + |A (1) , fin n (cid:105) , (2.9) |A (2) n (cid:105) = Z (2) |A (0) n (cid:105) + Z (1) |A (1) , fin n (cid:105) + |A (2) , fin n (cid:105) . (2.10)Here |A ( L ) n (cid:105) is a vector in colour space and admits the following expansion in the strongcoupling |A n (cid:105) = g n − s (cid:26) |A (0) n (cid:105) + α s π |A (1) n (cid:105) + (cid:16) α s π (cid:17) |A (1) n (cid:105) + O ( α s ) (cid:27) . (2.11)We remind that for the leading colour ¯ ttgg amplitudes the colour decomposition is givenin Eq. (2.2). The IR divergences of the amplitude are encoded in the Z factor (which is amatrix in colour space) and up to two-loop order. For the process we are considering, the Z factor is given by [55] Z = 1 + α s π (cid:18) Γ (cid:48) (cid:15) + Γ (cid:15) (cid:19) + (cid:16) α s π (cid:17) (cid:26) (Γ (cid:48) ) (cid:15) + Γ (cid:48) (cid:15) (cid:18) Γ − β (cid:19) + Γ (cid:15) (cid:0) Γ − β (cid:1) + Γ (cid:48) (cid:15) + Γ (cid:15) (2.12) − T F n h (cid:88) i =1 (cid:20) Γ (cid:48) (cid:18) (cid:15) ln µ m i + 14 (cid:15) (cid:20) ln µ m i + π (cid:21)(cid:19) + Γ (cid:15) ln µ m i (cid:21)(cid:27) + O ( α s ) . In this work we will consider the top quark as the only massive fermion, therefore n h = 1 and m = m t . The coefficients Γ n and Γ (cid:48) n are defined through the following expansion in α s , Γ = (cid:88) n ≥ Γ n (cid:18) α s π (cid:19) n +1 and Γ (cid:48) = − C i γ cusp ( α s ) . (2.13)For the ¯ ttgg amplitude, C i = C A = N c . The anomalous dimension matrix Γ for the ¯ ttgg amplitude can be obtained from Eq. (55) of [55]. We also refer the reader to Appendix Aof [52] for the definition of anomalous dimension factors required for the computation ofthe IR poles of the ¯ ttgg amplitude. There are a variety of different methods on the market for dealing with the spin states of amassive fermion. The common aim is to find a compact notation at the amplitude level andretain all information required to account for decays in the narrow width approximation.– 5 –ne major issue is that the helicity for a massive particle is not conserved and so we mustintroduce an arbitrary additional direction in order to define it. This additional extradirection can increase the algebraic complexity of an analytic computation.In this section we review one way in which massive spinors can be incorporated into thespinor-helicity methods [56] and then describe how to define a set of gauge invariant, on-shell sub-amplitudes that can be used to describe the full set of spin correlated narrow widthdecays. These on-shell sub-amplitudes can be computed using a rational parametrisation ofthe external kinematics. In our case, this rational paramtrisation was generated using themomentum twistor formalism [57]. The spinor-helicity method applied to massive fermionshas been well studied and we refer the reader to Refs. [58–60] for further details and otherapproaches.The first step is to introduce an arbitrary direction n which can be used to define amassless projection of the massive fermion momentum, p (cid:91),µ = p µ − m p.n n µ , (3.1)where p µ is the momentum of an on-shell massive fermion with p = m . The momenta p (cid:91) and n are both massless. The set of u and v spinors for the massive fermion can then beconstructed using the Weyl spinors for the massless momenta, u + ( p, m ) = ( /p + m ) | n (cid:105)(cid:104) p (cid:91) n (cid:105) , u − ( p, m ) = ( /p + m ) | n ][ p (cid:91) n ] , (3.2) v − ( p, m ) = ( /p − m ) | n (cid:105)(cid:104) p (cid:91) n (cid:105) , v + ( p, m ) = ( /p − m ) | n ][ p (cid:91) n ] . (3.3)The fact that the helicity state depends on the choice of reference vector means that positiveand negative helicities are no longer independent as they are in the massless case: u − ( p, m ) = (cid:104) p (cid:91) n (cid:105) m (cid:32) u + ( p, m ) (cid:12)(cid:12)(cid:12)(cid:12) p (cid:91) ↔ n (cid:33) . (3.4)For a t ¯ t pair this means that we only need to compute one helicity configuration with generalreference vectors, say ++ , and we will obtain enough information for all four possible spinconfigurations.We can then organise amplitudes for the ++ configuration into a basis of spinor struc-tures which parametrise the dependence on the two reference vectors n and n . Thereare four independent terms in this basis, matching up with the four configurations in the– 6 –omplete system. For the L -loop ¯ ttgg amplitude this is A ( L ) (1 +¯ t , + t , h , h ; n , n ) = m Φ h h (cid:104) (cid:91) n (cid:105)(cid:104) (cid:91) n (cid:105) (cid:32) (cid:104) n n (cid:105) A ( L ) , [1] (1 +¯ t , + t , h , h )+ (cid:104) n (cid:105)(cid:104) n (cid:105)(cid:104) (cid:105) A ( L ) , [2] (1 +¯ t , + t , h , h )+ s (cid:104) n (cid:105)(cid:104) n (cid:105)(cid:104) | | (cid:105) A ( L ) , [3] (1 +¯ t , + t , h , h )+ s (cid:104) n (cid:105)(cid:104) n (cid:105)(cid:104) | | (cid:105) A ( L ) , [4] (1 +¯ t , + t , h , h ) (cid:33) (3.5)where Φ is a phase factor depending on the helicities of the gluons. There are a few thingsto note about this expression:1. The sub-amplitudes A ( L ) , [ a ] are defined to be dimensionless and free of the spinorphases.2. A ( L ) , [ a ] depends only on two variables.3. The gluon phases for the ++ and + − configuration are taken to be Φ ++ = [34] (cid:104) (cid:105) , Φ + − = (cid:104) | | (cid:104) | | .
4. The basis of spin structures has been chosen to prefer (cid:104) n n (cid:105) over (cid:104) n (cid:105)(cid:104) n (cid:105) . Thesub-amplitudes appear to be simpler in this basis.The computation of the sub-amplitudes can then be obtained from four different eval-uations of the full amplitude with a rational kinematic configuration with four choices ofthe reference vectors. These evaluations can then be used to make a linear system which issolved to find the sub-amplitudes A ( L ) , [ a ] . Explicitly, A ( L ) (3 ,
3) = m Φ h h s (cid:104) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) | | A ( L ) , [4] , (3.6) A ( L ) (4 ,
4) = − m Φ h h s (cid:104) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) | | A ( L ) , [3] , (3.7) A ( L ) (3 ,
4) = m Φ h h (cid:104) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) (cid:91) (cid:105) A ( L ) , [1] , (3.8) A ( L ) (4 ,
3) = − m Φ h h (cid:104) (cid:105)(cid:104) (cid:91) (cid:105)(cid:104) (cid:91) (cid:105) ( A ( L ) , [1] + A ( L ) , [2] ) , (3.9)where we have dropped the particle labels on the A functions for simplicity. We start from a rational configuration with six massless particles generated via momentumtwistors [57, 61]. This configuration will depend on eight independent parameters. The– 7 –ix particles are then used to form a configuration with specific choices for the referencemomenta. We label the six massless momenta as q , . . . , q and the resulting t ¯ t system as p , p , p , p where, p = q + q , p = q + q , p = q , p = q , (3.10) q · q = q · q , (cid:104) q q (cid:105) = 0 , [ q q ] = 0 , (cid:104) q q (cid:105) = 0 , [ q q ] = 0 . (3.11)The ordering of the momenta q , . . . , q is important to find a simple rational solutionto the last four constraints. From the on-shell momenta p , . . . , p it is straightforwardto generate rational spinors for p (cid:91) , p (cid:91) using four different choices of reference vectors: ( n , n ) = { (3 , , (4 , , (3 , , (4 , } . The variables of the original six-point configuration,i.e. those of the q i momenta, are changed so that we can use conventional Mandelstaminvariants, s = ( p + p ) = 2 p · p , t = ( p + p ) = 2 p · p + m t . (3.12)This procedure would also work for high multiplicity amplitudes though a careful choiceof variables in the parametrisation may be necessary to obtain manageable algebraic com-plexity.The final results for the sub-amplitudes in Eq. (3.5) will be expressed using kinematicvariables x and y defined by − sm t = (1 − x ) x , tm t = y. (3.13)This choice rationalises the square root, (cid:113) s ( s − m t ) that appears in the master integral and amplitude computations. For the ¯ ttgg amplitudesthat do not involve elliptic sectors, the rationalisation of such a square root allows us toperform a computation within finite field arithmetic without the need to introduce addi-tional variables. It also important to make a rational change of variables when solving thedifferential equations of the master integrals. Following the procedure above allows us to directly evaluate the relevant colour orderedFeynman diagrams to obtain simple expressions at tree level: A (0) , [1] (1 +¯ t , + t , + , + ) = − − y ,A (0) , [2] (1 +¯ t , + t , + , + ) = A (0) , [3] (1 +¯ t , + t , + , + ) = A (0) , [4] (1 +¯ t , + t , + , + ) = 0 , (3.14)and A (0) , [1] (1 +¯ t , + t , + , − ) = A (0) , [3] (1 +¯ t , + t , + , − ) = ( x − y )(1 − xy )(1 − x ) (1 − y ) ,A (0) , [2] (1 +¯ t , + t , + , − ) = A (0) , [4] (1 +¯ t , + t , + , − ) = 0 . (3.15)– 8 –he procedure for generating loop level expressions does not change although the expres-sions require, as usual, tensor integral reduction to be performed. Our implementation ofthese steps is given in the subject of the following sections. Our final aim is to presentthe finite remainders of the sub-amplitudes after the subtraction of infrared and ultravioletpoles. The reduction of the helicity amplitudes, generated using Feynman diagrams, is performeddirectly 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 withfinite field arithmetic for problems with many external scales where the algebraic complexityis too high for conventional approaches. There are however many potential bottlenecks inamplitude reduction and it has been shown to be important to control as many steps inthe computation as possible. Amplitudes requiring dimensional regularisation can causeunnecessary large intermediate steps since information that vanishes in four dimensions isretained. Unitarity cut based approaches [62] and specially designed tensor projectors [63–65] have been used successfully in recent high multiplicity loop amplitude computations inmassless theories.In our setup the diagram generation, colour ordering and spinor-helicity algebra areperformed with the help of
Qgraf [66],
Form [67, 68],
Spinney [69] and
Mathematica scripts. During this phase we identify a set of topologies that are independent at theintegrand level. The numerators of these topologies are then separated into loop momentumdependent structures and coefficients depending on the external kinematics. The latter arethen evaluated using the rational parametrisation described in the previous section. Up tothis point we follow a strategy described in, for example, Ref. [70]. The only difference isthat the internal masses must also be tracked.The result of these steps is not yet suitable for processing with integration-by-partsidentities since the loop dependent numerators must be transformed into a basis of inde-pendent scalar products thereby upgrading the topologies to complete families of Feynmanintegrals. In the case of four-particle amplitudes an additional integration over the spuriousspace in the loop momenta must also be performed. We employ a transverse integrationapproach similar to the one outlined in Ref. [71].For clarity, we present some additional details of these steps that are specific to ourproblem. An interested reader may wish to refer to Refs. [72–74] for a more completediscussion and other applications. All steps described from this point to the reduction ontoa basis of special functions have been constructed using
FiniteFlow graphs. After colourordering and helicity amplitude processing the amplitude takes the following form, A ( L ) ,h ( { p } ) = (cid:90) L (cid:89) j =1 d d k j iπ d/ e − (cid:15)γ E (cid:88) T N hT ( d s , { k } , { p } ) (cid:81) α ∈ T D α ( { k } , { p } ) , (4.1)where d = 4 − (cid:15) , d s = g µµ and k i are the spacetime dimension, spin dimension and loopmomenta, respectively. T is a set of distinct diagram topologies contributing to the colour-– 9 – T T T T T T T T
123 41243 T T Figure 4 : Maximal topologies for the N c part of planar two-loop ¯ ttgg amplitude. Black-solid lines represent top quarks, red-dashed lines represent gluons.stripped helicity amplitude A ( L ) ,h ( { p } ) . A diagram topology may include contributionsfrom more than one Feynman diagram. To parameterise the loop integrand, we decomposethe d -dimensional loop momenta k i into a four-dimensional part ¯ k i and an extra-dimensionalpart ˜ k i , k µi = ¯ k µi + ˜ k µi . (4.2)Since the external momenta are in four dimensions, the extra dimensional part of the loopmomenta only appears in the numerator function as µ ij = − ˜ k i · ˜ k j . After substitutingmassive spinors in terms of massless ones according to Eqs. (3.2) - (3.3) and performingt’Hooft algebra, the numerator function contains the following loop-momentum dependentterms k i · k j , ¯ k i · p j , µ ij , (cid:104) p a | ¯ k i | p b ] , (cid:104) p a | ¯ k i ¯ k j | p b (cid:105) , [ p a | ¯ k i ¯ k j | p b ] , (4.3)where momenta p a and p b are massless and p j can either be massive or massless. As indi-cated above, to apply IBP reduction we must first express these loop momentum structuresin a basis of the propagators D α and a set of irreducible scalar products (ISPs) which definea family of Feynman integrals.We make sure that we introduce the minimal number of integral families by mappingall distinct diagram topologies, T , to a set of master topologies which have the maximumnumber of propagators. The master topologies for the N c part of the two-loop planar ¯ ttgg 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 parts ¯ k µ (cid:107) ,i and spurious parts ¯ k µ ⊥ ,i , ¯ k µi = ¯ k µ (cid:107) ,i + ¯ k µ ⊥ ,i , (4.4)– 10 –nd we further expand them into a physical ( v µ ) and a spurious ( w µ ) spanning bases ¯ k µ (cid:107) ,i = d (cid:107) (cid:88) j =1 a ij v µj , ¯ k µ ⊥ ,i = d ⊥ (cid:88) j =1 b ij w µj . (4.5) d (cid:107) = min(4 , n ext ) is the dimension of the d physical space, n ext is the number of indepen-dent momenta of a given topology and d ⊥ = 4 − d (cid:107) is the dimension of the d spuriousspace. 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 andphysical ISPs, a ij = a ij ( D α , ISPs) , while the coefficients of the spurious spanning basis arefunctions of spurious ISPs and extra dimensional ISPs, b ij = b ij ( k i · w j , µ ij ) . a ij and b ij aredetermined by solving a linear system of equations obtained by contracting Eq. (4.5) with v i and w i . For the ¯ ttgg amplitudes, all the maximal topologies have four external momenta,hence n ext = 3 , d (cid:107) = 3 and d ⊥ = 1 , and for instance, we can choose the following set ofspanning vectors v µ = { p µ , p µ , p µ } , w µ = { ω µ } . (4.6)Explicitly, we used: ω µ = 12 (cid:0) (cid:104) | | | σ µ | − (cid:104) | σ µ | | | (cid:1) , (4.7)though any vector satisfying ω · p = ω · p = ω · p = 0 will suffice. Numerator termswith odd powers of k i · w j vanish under loop integration, while the terms with even powerscan be cast into linear combination of propagator denominators D α and physical ISPs. Forexample, k i · w k k j · w k = w k d ⊥ − (cid:15) (cid:0) k i · k j − k (cid:107) ,i · k (cid:107) ,j (cid:1) , (4.8)where k i · k j can be directly written in terms of D α while k (cid:107) ,i · k (cid:107) ,j can be expressed in termsof D α and physical ISPs according to Eq. (4.5). Similarly for the extra dimensional ISP, µ ij = − (cid:15)d ⊥ − (cid:15) (cid:0) k i · k j − k (cid:107) ,i · k (cid:107) ,j (cid:1) . (4.9)After the replacements have been made the amplitude is recast into a form suitable for IBPreduction, A ( L ) ,h ( { p } ) = (cid:88) T (cid:88) i c hT,i ( (cid:15), { p } ) G T,i ( (cid:15), { p } ) . (4.10)Note that there may be overlap in sets of Feynman integrals G T,i ( { k } , { p } ) appearing ineach numerator topology T . This is in contrast to the integrand reduction approach takenin [70] where there is no overlap in G T,i ( { k } , { p } ) between numerator topologies. Combiningcontributions from all numerator topologies we have A ( L ) ,h ( { p } ) = (cid:88) k ˜ c hk ( (cid:15), { p } ) ˜ G k ( (cid:15), { p } ) , (4.11)where ˜ G 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– 11 –ith 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 am-plitude, A ( L ) ,h mren ( { p } ) = A ( L ) ,h ( { p } ) + A ( L ) ,h mct ( { p } ) , (4.12) = (cid:88) k c mren ,hk ( (cid:15), { p } ) ˜ G k ( (cid:15), { p } ) . (4.13)The mass renormalisation counterterms appearing in Eqs. (2.5) and (2.6) have to be writtenin 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 integralsvia IBP identities. The IBP relations are generated in Mathematica with the help of
LiteRed [75], and solved numerically over finite fields within the
FiniteFlow frame-work [38] using the Laporta algorithm [76]. The mass renormalised helicity amplitude isnow written in a basis of master integrals MI k , A ( L ) ,h mren ( { p } ) = (cid:88) k c IBP ,hk ( (cid:15), { p } ) MI k ( (cid:15), { p } ) . (4.14)Let us note that not all of the maximal topologies that appear in an amplitude are processedin the IBP system, as some of them can be mapped into the other maximal topologies. Forexample, maximal topologies T , . . . , T in Fig. 4 can be mapped to topologies T , . . . , T by choosing an appropriate set of physical ISPs (auxiliary propagators) for each of thosetopologies. The master integral representation can be used to check Ward identities onthe mass-renormalised amplitude by replacing gluon polarisation with its momentum inthe numerator construction and reconstructing Eq. (4.14) over a prime field to check if itindeed vanishes.If the analytic solutions of the master integrals appearing in Eq. (4.14) for the pla-nar ¯ ttgg amplitudes are available, we can proceed further by writing the master integralsas a linear combination of special function monomials m ( f ( { p } )) and perform Laurentexpansions, A ( L ) ,h mren ( { p } ) = (cid:88) k (cid:88) l = n ( L ) (cid:15) l c exp ,hkl ( { p } ) m k ( f ( { p } )) + O ( (cid:15) ) , (4.15)where n ( L ) is the power of the deepest pole that can appear in the L -loop amplitude. Forone- and two-loop cases, n (1) = − and n (2) = − . To obtain the simplest representation ofthe amplitude we further subtract the UV and IR poles from the mass-renormalised ampli-tude. 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 planar ¯ ttgg one- and two-loop amplitudes, in the same special function mono-mial basis as in Eq. (4.15), P ( L ) ,h UV ( { p } ) = (cid:88) k (cid:88) l = n ( L ) (cid:15) l c UV ,hkl ( { p } ) m k ( f ( { p } )) + O ( (cid:15) ) , (4.16) P ( L ) ,h IR ( { p } ) = (cid:88) k (cid:88) l = n ( L ) (cid:15) l c IR ,hkl ( { p } ) m k ( f ( { p } )) + O ( (cid:15) ) , (4.17)– 12 –nd subtract them from the mass-renormalised amplitude to obtain the finite remainder, F ( L ) ,h ( { p } ) = A ( L ) ,h mren ( { p } ) + P ( L ) ,h UV ( { p } ) − P ( L ) ,h IR ( { p } ) (4.18) = (cid:88) k c F ,hk ( { p } ) m k ( f ( { p } )) + O ( (cid:15) ) . (4.19)The coefficients c F ,hk ( { p } ) appearing in (4.19), however, are not all independent. We exploitthis fact in order to simplify both the result and the reconstruction of its analytic expressionas follows. First, we sort all the coefficients by their complexity, which is estimated fromtheir total degree. The total degree, in turn, can be quickly determined via a univariatefit, as explained in [37]. We then find vanishing linear combinations of these coefficients bysolving the linear fit problem (cid:88) k y k c F ,hk ( { p } ) = 0 (4.20)with respect to the unknowns y k , which are rational numbers. This allows us to find linearrelations between the coefficients, which rewrite the more complex ones in terms of simplerones. After applying the linear relations between coefficients of special function monomials,we arrive at F ( L ) ,h ( { p } ) = (cid:88) k ¯ c F ,hk ( { p } ) ¯ m k ( f ( { p } )) + O ( (cid:15) ) , (4.21)where ¯ c F ,hk are the independent coefficients of the new special function monomials ¯ m k ( f ) ,which are linear combinations of the monomials appearing in Eq. (4.19). Therefore, func-tional reconstruction only needs to be applied to the independent coefficients ¯ c F ,hk . Thisyields a significantly simpler result than the one in Eq. (4.19), and also reduces the numberof evaluations needed for its reconstruction.Having set up the numerical algorithm suitable for computation over finite fields, start-ing from the Feynman diagram numerators, to evaluate ¯ c F ,hk in the L -loop finite remainder F ( L ) ,hn ( { p } ) using the FiniteFlow framework, the analytic forms of ¯ c F ,hk are obtained byusing the multivariate reconstruction method described in [37], after performing severalnumerical evaluations. To derive an analytic representation of the helicity amplitudes, particularly when we areinterested in computing the finite remainder, analytic expressions of master integrals ap-pearing in the one- and two-loop amplitudes are required. The solutions of all one- andtwo-loop master integrals appearing in the A (2) , , A (2) ,N l , A (2) ,N l , A (2) ,N l N h and A (2) ,N h amplitudes in Eq. (2.4) can be expressed in terms of multiple polylogarithms (MPLs) [81–83] and have been completed recently [27, 77]. For the two-loop amplitude A (2) ,N h involvinga single top-quark closed-loop, the master integrals also involve elliptic generalisations ofMPLs [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 am-plitudes: – 13 – a ) ( b ) ( c ) Figure 5 : Master integral topologies with maximal number of propagators (7) that appearin the leading colour 2-loop amplitude for ¯ ttgg . Black-solid lines represent massive particles,red-dashed lines represent massless particles. Figure 6 : New master integral topologies appearing in the amplitude with a closed top-quark loop A (2) ,N h . Black-solid lines represent massive particles, red-dashed lines representmassless particles.• A (2) , , A (2) ,N l , A (2) ,N l : There are 36 master integrals appearing in the A (2) , ampli-tude. We have used the results in [77] for the relevant master integrals. The basisof master integrals involves the topologies with 7 propagators (double-box) shown inFig. 5 (topology a and b ). The master integrals required for the computation of A (2) ,N l and A (2) ,N l can be obtained from a subset of the master integral basis of A (2) , .• A (2) ,N h : There are 55 master integrals in the amplitude with a single top-quark closedloop. 44 of them belong to the topbox family, i.e. the planar double-box integralfamily contributing to A (2) ,N h (topology c in Fig. 5), computed in [18]. Results forthe 9 master integrals that are not part of the topbox family are obtained from [77].Analytic expressions for the remaining 2 master integrals (shown in Fig. 6), to the bestof our knowledge, are not available in the literature. In this work we have computedthose two integrals using the method of differential equations. The details of thecomputation will be covered in the next subsection.• A (2) ,N l N h and A (2) ,N h : Both amplitudes are of one-loop squared type, and the masterintegrals needed are the one-loop bubbles and triangle up to weight 4, which arealready available from the one-loop amplitudes computation up to O ( (cid:15) ) .– 14 – ( p )¯ t ( p ) t ( p ) g ( p ) ttt Figure 7 : Feynman diagram contributing to A (2) ,N h that leads to the new master integraltopologies shown in Fig. 6. The solid lines denoted the top quark whereas spiral lines denotethe gluons. All external momenta are on-shell. The new master integral topologies that appear in A (2) ,N h arise from the Feynman diagramshown in Fig. 7. The integral family for this penta-triangle topology is given by I ν ν ν ν ν ν ν ν ν = e (cid:15)γ E (cid:0) m t (cid:1) ν − d (cid:90) d d k iπ d/ d d k iπ d/ D − ν D − ν D ν D ν D ν D ν D ν D ν D ν , (5.1)where γ E denotes the Euler-Mascheroni constant, ν = (cid:80) j =1 ν j and d = 4 − (cid:15) . The inversepropagators D i are given by D = − k , D = − ( k − p ) , D = − ( k + p + p ) + m t ,D = − ( k + p ) , D = − k + m t , D = − ( k − p ) + m t ,D = − ( k + k ) + m t , D = − ( k + p ) , D = − ( k + p + p ) . (5.2)We choose a basis of master integrals for which the irreducible scalar products correspondingto D and D are absent. Therefore we may label these master integrals by I ν ν ν ν ν ν ν .Using IBP identities, we obtain 22 master integrals, { M i } 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 M = I , M = I , M = I ,M = I , M = I , M = I ,M = I , M = I , M = I ,M = I , M = I , M = I ,M = I , M = I , M = I ,M = I , M = I , M = I ,M = I , M = I , M = I ,M = I . (5.3)– 15 –
23 4 M M M M M
532 41 M M M M M
102 3 4 1 M
11 3 412 M
12 32 14 M
13 3 412 M
14 3 412 M M
16 4123 M
17 32 14 M
18 32 14 M
19 3 42 1 M M
21 4123 3 2 14 M Figure 8 : Master integral basis for the pentatriangle topology shown in Fig. 7. Black-solidlines represent massive particles, red-dashed lines represent massless particles.A set of 15 integrals were identified as master integrals from the topbox family and werealready computed in [18]. They are M = I , M = I , M = I , M = I ,M = I , M = I , M = I , M = I ,M = I , M = I , M = I , M = I ,M = I , M = I , M = I , (5.4)where I ν ν ν ν ν ν ν are the precanonical master integrals defined in Eq. (8) of [18]. Asecond set of 5 master integrals is available from [77], and they are identified as M = T , M = T , M = T , M = T , M = T , (5.5)– 16 –here T i are master integrals defined in Sec. 6.2 of [77]. M and M are the missingintegrals. Note that these two integrals also appear in a planar double-box family thatcontributes to the subleading colour ¯ ttgg amplitude without closed fermion loops (topology P in Figure 1 of [78]).Apart from the x and y defined in Eq. (3.13), it can be useful to use another coordinatesystem ˜ x and ˜ y defined as − sm t = (1 + ˜ x ) ˜ x (1 − ˜ x ) , y = 1 − ˜ y , (5.6)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 d (cid:126) J = (cid:15)A (cid:126) J , A = A x dx + A y dy, (5.7)where A does not depend on (cid:15) 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 familyis chosen as J = − (cid:15) yM , J = (cid:15) (1 − y ) M , J = − (cid:15) (1 − (cid:15) ) (3 − (cid:15) )(1 − y ) M − (cid:15) (1 − (cid:15) )(2 − (cid:15) )(1 − (cid:15) )(1 − y ) M + 12 (cid:15) (1 − y )(1 + y ) M , (5.8) J = − (cid:15) (1 − (cid:15) )(1 − y ) M − (cid:15) (1 − y )(1 + x ) x M , J = (cid:15) (1 − y )(1 − x )(1 + x )2 x M , J = − (cid:15) (1 − y ) M , J = − (cid:15) (1 − x ) (1 − y ) x M . For the set of master integrals belonging to the topbox topology we follow the choice givenin Sec. 6 of [18]. Here the form of the differential equation is slightly relaxed and thedifferential equation for the basis J has the form d (cid:126) J = A (0) x,y + (cid:15)A (1) x,y (cid:126) J , (5.9)where A (0) is strictly lower triangular and A (0) and A (1) are independent of (cid:15) . A ( i ) x,y is rationalin the kinematic variables ( x, y ) , periods of elliptic curves and their derivatives [21]. Byconsidering them as independent variables, we derived the system of differential equationsefficiently with the finite field sampling method using the FiniteFlow package [38]. Wenote that in the computation of the A (2) ,N h finite remainder, the coefficients of the special– 17 –unction monomials, ¯ c F ,hk in Eq. (4.21), are also rational functions in ( x, y ) , periods of ellipticcurves and their derivatives. The solution for master integrals in terms of iterated integralsare presented up to 4 orders in (cid:15) : J i = (cid:88) k =0 (cid:15) k J ( k ) i + O ( (cid:15) ) . (5.10) Now that we have identified the necessary master integrals, let us discuss the mathematicalobjects 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. Letus therefore start by reviewing the definition of Chen’s iterated integrals [79]. Considera path γ on an n -dimensional manifold M with starting point x i = γ (0) and end point x f = γ (1) , γ : [0 , → M. Let us also consider a set of differential 1-forms { ω i } and their pullbacks to the interval [0 , which we denote by f j ( λ ) dλ = γ ∗ ω j . (5.11)Then we define the k -fold iterated integral over ω , ..., ω k along γ as I γ ( ω , ..., ω k ; λ ) = (cid:90) λ dλ f ( λ ) (cid:90) λ dλ f ( λ ) ... (cid:90) λ k − dλ k f k ( λ k ) (5.12) = (cid:90) λ dλ f ( λ ) I γ ( ω , ..., ω k ; λ ) . (5.13)From the above definition it is clear that the derivative of an iterated integral is given by ddλ I γ ( ω , ω , ..., ω k ; λ ) = f ( λ ) I γ ( ω , ..., ω k ; λ ) . It can be shown that the iterated integrals defined in Eq. (5.12) have various interestingproperties [80]. The most useful property for the purpose of this paper is concerning thedecomposition of the path γ . Let α, β : [0 , → M be two paths such that β (0) = α (1) andlet γ = αβ be the path obtained by concatenating α and β . Then we can decompose I γ ( ω , ..., ω k ; λ ) = n (cid:88) i =0 I β ( ω , ..., ω i ; λ ) I α ( ω i +1 , ..., ω n ; λ ) , (5.14)where we define the -fold integrals as I γ (; λ ) = 1 . (5.15)– 18 –ow that we are familiar with the notion of iterated integrals, let us consider the well-understood multiple polylogarithms (MPLs). MPLs are a special class of iterated integralswhere the -forms ω i are such that γ ∗ ω j = dλλ − c , (5.16)for some c ∈ C . Then we have G ( z , ..., z k ; y ) = (cid:90) y dy y − z G ( z , ..., z k ; y ) . (5.17)Note that the recursion defined in (5.17) diverges whenever z k = 0 . We can regularise thisdivergence by setting G (0 , . . . , (cid:124) (cid:123)(cid:122) (cid:125) k times ; y ) = 1 k ! (ln y ) k . (5.18)The class of MPLs is very well understood and many tools for their algebraic manipu-lation and their numerical evaluation are available [84–86], it is however not always possibleto express all master integrals in terms of these functions. Especially in computations in-volving massive internal particles we often encounter elliptic sectors which transcend thespace of MPLs. In the following, we give a brief overview of the elliptic sectors we encounterin this work as well as the computational peculiarities that come along with them. In this section we review some basic properties of elliptic curves as well as the elliptic curvesspecific to the topbox topology, which were previously analysed in [18]. For more thoroughreviews of elliptic curves and their properties, we refer the reader to one of the variousanalyses in the literature [18, 87, 88]. Let us start by considering the generic quartic caseof an elliptic curve. An elliptic curve E can be described by the equation E : w − ( z − z ) ( z − z ) ( z − z ) ( z − z ) = 0 , (5.19)where the z i ∈ C . The z i define the properties of the elliptic curve and are generallyfunctions of the kinematic variables x = ( x , ..., x n ) : z j = z j ( x ) , j ∈ { , , , } . (5.20)Alternatively, an elliptic curve can be described in terms of two elliptic periods ψ i . In orderto define these elliptic periods, let us introduce the auxiliary variables Z i as Z = ( z − z ) ( z − z ) , Z = ( z − z ) ( z − z ) , Z = ( z − z ) ( z − z ) . (5.21)The Z i satisfy the relation Z + Z = Z (5.22)– 19 – a) The sunrise topology (b) The topbox topology (c) The bubblebox topology Figure 9 : Topologies in the topbox family that are associated with the elliptic curves E ( a ) , E ( b ) and E ( c ) . Black-solid lines represent massive particles, red-dashed lines representmassless particles.and can be used to define the modulus k and the complementary modulus ¯ k of the ellipticcurve E as k = Z Z , ¯ k = 1 − k = Z Z . (5.23)Then we can choose the two elliptic periods ψ i associated to the elliptic curve E as ψ = 4 K ( k ) Z , ψ = 4 iK (cid:0) ¯ k (cid:1) Z , (5.24)where K denotes the complete elliptic integral of the first kind K ( λ ) = (cid:90) dt (cid:112) (1 − t )(1 − λt ) . (5.25)In our results, the dependence of iterated integrals on elliptic curves enters throughthe appearance of elliptic periods in the integration kernels Eq. (5.11). Let us now take acloser look at the different elliptic topologies and the corresponding elliptic curves relevantfor this publication. The topbox diagram has three elliptic sub-sectors corresponding tothree different elliptic curves. They can be obtained from the maximal cuts in the Baikovrepresentation [89]. 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 probablythe most well-known of the three. The corresponding z i in Eq. (5.19) are given as z ( a )1 = t − m µ , z ( a )2 = − m − m √ tµ , z ( a )3 = − m + 2 m √ tµ , z ( a )4 = tµ . (5.26)Iterated integrals involving only this topology can be cast in terms of iterated integrals ofmodular forms, and are well-suited for numerical evaluation [19, 23, 24, 90–93].The elliptic curve E ( b ) is associated to the topbox sector itself (Fig. 9b) and the corre-sponding z i are given by z ( b )1 = t − m µ , z ( b )2 = − m − m (cid:113) t + ( m − t ) s µ , z ( b )3 = − m + 2 m (cid:113) t + ( m − t ) s µ ,z ( b )4 = tµ . (5.27)– 20 –he elliptic curve E ( c ) is associated to the remaining elliptic sub-sector which we dubthe bubblebox sector (Fig. 9c). The corresponding quartic is defined by z ( c )1 = t − m µ , z ( c )2 = 1 µ (cid:32) − m ( s + 4 t )( s − m ) − m − s (cid:114) sm (cid:16) st + ( m − t ) (cid:17)(cid:33) ,z ( c )3 = 1 µ (cid:32) − m ( s + 4 t )( s − m ) + 24 m − s (cid:114) sm (cid:16) st + ( m − t ) (cid:17)(cid:33) , z ( c )4 = tµ . A more thorough analysis of the different sub-sectors and the corresponding elliptic curvescan be found in [18].
Now that we have identified the various elliptic curves entering the computation, let usturn to the computation of the master integrals J and J , which are not known from theliterature. As we have previously explained, the missing master integrals can be computedorder by order in (cid:15) from the differential equation given in Eq. (5.9). We compute the masterintegrals iteratively in the expansion parameter (cid:15) by integrating the differential equationfrom a base point ( x , y ) = (0 , (corresponding to s = ∞ and t = m ). The integralsthat need to be evaluated receive contributions from the elliptic topbox sub-sectors andhence contain elliptic iterated integrals themselves. The boundary values of the integralsat different orders in (cid:15) are inferred from the regular point x, y = (1 , , where these twomaster integrals vanish.The integral J is associated with elliptic curve a , and the first five orders in (cid:15) for J read J (0)21 = 0 , J (1)21 = 0 , J (2)21 = 0 , J (3)21 = 0 , J (4)21 = − π
60 + (cid:0) G (0 , y ) − G (1 , y ) (cid:1) ζ + G (0 , , , , y ) − G (1 , , , , y )+ π I γ ( g , f ; λ ) + 118 I γ ( g , f , η ( a )0 , f ; λ ) . (5.28)We have adopted the notation for the integration kernels appearing in the topbox familywhich was introduced in [18]. The integration kernels in Eq. (5.28) are given by g = dy y + 3 y (1 − y ) , (5.29) f = 3 ψ ( a )1 π dy , (5.30) η ( a )0 = − π dy ( ψ ( a )1 ) ( − y )( − y ) y . (5.31)– 21 –e find that J contains the elliptic sector 93 from [18] in one of its sub-topologiesand is hence associated with the elliptic curve b . The first four orders in (cid:15) for J can bewritten entirely in terms of MPLs and read J (0)22 = 0 , J (1)22 = 0 , J (2)22 = − G (0 , , x ) , J (3)22 = G (1 , y ) G (0 , , x ) + 3 G (0 , − , , x ) − G (0 , , , x ) + π G (0 , x ) + 92 ζ . (5.32)The explicit dependence on the elliptic curve associated with the topbox sector enters forthe first time at O ( (cid:15) ) . The corresponding term involves 28 new integration kernels thathave not appeared in the computations of the other master integrals and, accordingly, yieldsa rather large expression. We therefore refrain from showing it here explicitly and refer theinterested reader to the ancillary files.Note that the differential -forms appearing in the iterated integrals are generally notpath independent [80]. Nevertheless it is clear from a physical point of view that the masterintegrals do not depend on the particular choice of path chosen for the numerical evaluationof these iterated integrals. We have verified this path independence numerically using themethods explained in the following section. In the previous sections, we discussed the analytic expressions of all the master integralsrelated to our process. In this section, we discuss how to numerically compute these inte-grals. The integrals containing only MPLs can be evaluated to a high precision, for example,using
Ginac [94]. The numerical evaluation of iterated integrals associated with ellipticcurves on the other hand is quite challenging. We explain how to compute these iteratedintegrals 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 theboundary point is relatively straightforward, as was also explained in [18]. For numericallyevaluating these integrals in the physical region we need to analytically continue the asso-ciated elliptic curves across all the branch cuts appropriately. The analytic continuationof integrals having one parameter dependence has already been discussed in the literature[22, 24, 93]. However, analytically continuing all the integrals in our case, which includesthe topbox integral associated to three different elliptic curves, is much more involved. Withthe objective of not digressing too much from the main goal of this paper, we discuss hereonly the ingredients essential for this computation.The physical region for our case is governed by the following equations ( m t = 1 ): s ≥ , G ( p , p , p ) ≥ ⇒ s ( − t − st + 2 t − ≥ , – 22 –here G is the Gram-determinant. In the coordinates x and y , these equations take theform y = t,x = 12 (2 − √ − s √− s − s ) . Fig. 10 shows the physical region in s, t coordinates as well as x, y coordinates. ��� ��� ��� ��� - � - � - � (a) Physical region in s t coordinates - ��� - ��� - ��� ��� ��� - � - � - ��� (b) Physical region in x y coordinates Figure 10 : The physical region for gg → t ¯ t .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 overthe whole path. The path segments need to be chosen while taking into consideration bothspurious as well as physical types of singularities in the system.The Euclidean region is defined as the region where the iterated integrals have onlyreal contributions. The numerical evaluation of the iterated integrals that appear in ourresults for the master integrals is done for the Euclidean region as follows. We split theintegration path from the boundary point (0 , to a generic kinematic point ( x, y ) in theEuclidean region into two pieces and use the path decomposition formula Eq. (5.14) toevaluate the integrals. It is better to use the coordinate system in Eq. 5.6 (˜ x, ˜ y ) here. ˜ x rationalizes all the square roots associated with the non-elliptic kernels and ˜ y is used tobring the boundary point of integration to (0 , . We choose the following paths: α , from (0 , to (˜ x, and β , from (˜ x, to (˜ x, ˜ y ) . For all the master integrals in our case, theintegration 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 integrationkernels around y = 1 , assuming ˜ y to be small. For integrating to a point y not close to ,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 continuethe 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 performthe integration. The choice of path is controlled by the radius of convergence of the series,– 23 –hich in turn depends on the singularities present in the kernels. It is important to chooseproperly the number of path segments and their sizes, so that the series solution convergesproperly on each of these segments and also the computation does not become too slowdue to the presence of too many segments. Generally, the segments should not be largerthan half the radius of convergence of the series expansion [20]. We again use the pathdecomposition formula to patch together the contribution from the different segments andcompute the iterated integral over the whole path. - ��� - ��� - ��� ��� ��� - � - � - ��� Figure 11 : Integrating to a point P in the physical region.For all the sunrise type kernels, there is a singularity while crossing the line y = 0 ( ˜ y = 1 ). For the elliptic kernels having both x and y dependence, there is a singularity onthe 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 theiterated integrals first along y and then along x . The kernels are no longer only MPLs oneither of these paths but evaluate to iterated integrals over elliptic curves.For all the integrals from the topbox family, we compute analytically in the physicalregion all the iterated integrals associated with the curves a and b , along with all theintegrals depending on MPLs only. The periods of elliptic curve c are discontinuous on theline x = 0 , unlike periods a and b . We reserve for a future publication a detailed explanationof the analytic continuation of the iterated integrals of the topbox family, which depend onmultiple elliptic curves, together with the study of the analytic continuation of the integralswhich depend also on the curve c .Numerical integration to other regions in the phase space, not discussed above, canalso be performed in a similar way. For example, integrating to the point ( − . , . inFig. 11 does not require us to cross any physical singularities if we integrate all the iteratedintegrals from the base point (0 , . This can also be performed in the same way using thepath decomposition formula recursively and does not involve a lot of complications. Anotherway to perform the integration in different phase space points is to start integrating theiterated integrals from a point already in that region. While this avoids the problem ofanalytic continuation it requires us to compute the boundary conditions at the new basepoints. – 24 – .1 Functional Relations Among Iterated Integrals In the previous sections, we have presented how we compute the planar master integralscontributing to the gg → t ¯ t scattering process in terms of MPLs and iterated integralsover elliptic curves. While this allows us to evaluate the t ¯ tgg helicity amplitudes, we findthat there are potential redundancies in the functional basis of iterated integrals that mightmake the expression appear complex. Indeed, after integrating the differential equation, wefound that the simple pole vanished numerically as expected but this cancellation was notreflected in the analytic expression. In order to achieve an explicit cancellation of all poleterms, some functional relations among the iterated integrals had to be applied. In thissection we comment on functional relations among iterated integrals and the cancellationof the (cid:15) pole. For this purpose, we adopt again the notation for the integration kernelsintroduced in [18].The functional relations we will discuss in this section are related to the appearance ofderivatives of elliptic periods in the integration kernels. Take for example the kernel a ( b )3 , = ψ ( b )1 xπ ( x − dy − ( ∂ y ψ ( b )1 ) x ( y − (cid:0) x y + 3 x − xy − x + y + 3 (cid:1) π ( x − ( x + 1) (3 x − xy − x + 3) dx (6.1) + ( ∂ y ψ ( b )1 ) x ( y − π ( x − dy − ψ ( b )1 ( y − (cid:0) x − x y + x − x y − xy + x + 3 (cid:1) π ( x − ( x + 1) (3 x − xy − x + 3) dx = (cid:32) ( ∂ x ψ ( b )1 ) x ( y − π ( x − − ψ ( b )1 ( x + 1)( y − π ( x − (cid:33) dx (6.2) + (cid:32) ( ∂ y ψ ( b )1 ) x ( y − π ( x − + ψ ( b )1 xπ ( x − (cid:33) dy = d (cid:18) ψ ( b )1 x ( − y ) π ( − x ) (cid:19) . (6.3)We see that the integration kernel a ( b )3 , corresponds to a total derivative. Because theprimitive of a ( b )3 , vanishes at the lower integration boundary, we have I ( a ( b )3 , ) = ψ ( b )1 x ( − y ) π ( − x ) . (6.4)We can furthermore use integration by parts to derive other relations for iterated inte-grals involving the kernel a ( b )3 , . In particular, we have I ( a ( b )3 , , f, . . . ) = (cid:90) a ( b )3 , I ( f, . . . ) (6.5) = (cid:90) d (cid:18) ψ ( b )1 x ( − y ) π ( − x ) (cid:19) · I ( f, . . . ) (6.6) = (cid:20) ψ ( b )1 x ( − y ) π ( − x ) I ( f, . . . ) (cid:21) ( x,y )(0 , − I (cid:18) ψ ( b )1 x ( − y ) π ( − x ) f, . . . (cid:19) , (6.7)– 25 –here f is an arbitrary integration kernel. Note that the first factor in the integral inEq. (6.6) denotes a total differential and not the integration measure. If we can recast theproduct ψ ( b )1 x ( − y ) π ( − x ) f (6.8)as a linear combination of integration kernels, we have successfully identified a relationamong iterated integrals I which is closed under the minimal set of integration kernels weare using. We find the following equalities ψ ( b )1 x ( − y ) π ( − x ) η ( b )0 = 23 η ( b )1 , , (6.9) ψ ( b )1 x ( − y ) π ( − x ) η ( a,b )2 = 172 η ( a )3 , − f , (6.10) ψ ( b )1 x ( − y ) π ( − x ) η ( b )1 , = 32 η ( b )2 , − ω + 18 ω . (6.11)Plugging the above linear combinations into Eq. (6.7), we obtain three functional relationswhich we can use to simplify the pole terms. After doing so we are left with iteratedintegrals involving only rational kernels. These integrals evaluate to MPLs and cancel theones 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 basisof 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, wefound that there were 12025 monomials which, applying relations involving a ( b )3 , as describedabove, was reduced to 11791 monomials. We notice that there are much fewer monomialbasis elements appearing in the two-loop finite remainders F (2) ,N h . For example, for oneof the subamplitudes with helicities + + + − , only 3586 monomial basis elements appear inthe finite remainder with non-zero coefficients. Applying the relations discussed above wefind that this number is reduced to 3158. Similar situations are also observed for the othersubamplitudes and the + + ++ helicity configuration.Note also that there are more kernels containing derivatives of elliptic periods. It ispossible to write some of these kernels as a sum of a total derivative which captures all ap-pearing derivatives of elliptic periods and a new kernel which does not contain the derivativeof an elliptic period. By doing so, it might be possible to find further relations. Because ofthe large number of iterated integrals contributing to the scattering amplitude and the needto introduce new integration kernels, from our experience with a ( b )3 , kernel, finding more re-lations between iterated integrals is necessary to simplify the finite remainder. We thereforereserve a thorough analysis of the possible functional relations for another publication.– 26 – Numerical Results
For benchmarking purposes, we present numerical results at a phase space point. We use aEuclidean phase space point that is considered in [18] p = ( − . , . , . i, . ,p = (75 . , . , − . i, − . ,p = (115 . , − . , . i, − . ,p = ( − . , . , − . i, . ,m t = 3 , (7.1)with the following choice of reference vectors η = (74 . , . , − . i, − . ,η = ( − . , . , . i, . . The rational representation of the momenta listed above is provided in the ancillary files.The set of momenta given in Eq. (7.1) corresponds to the following Mandelstam invariants s = − , t = 9011 , and x = 7120 , y = 1011 , ˜ x = 115 , ˜ y = 111 . We present results for the mass renormalised two-loop helicity amplitude for variousfermion loop contributions specified in Eq. (2.4), normalised to the tree-level amplitude ˆ A (2) ,fλ λ λ λ = A (2) ,f (1 λ ¯ t , λ t , λ g , λ g ) A (0) (1 λ ¯ t , λ t , λ g , λ g ) , (7.2)with helicities λ i and f = 1 , N l , N h , N l , N l N h , N h . Numerical results at the kinematic pointgiven in Eq. (7.1) are displayed in Table 1. We use Ginac to evaluate MPLs that appear inthe amplitude, while the numerical evaluation of iterated integrals involving elliptic kernelsare 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 elementobtained from the helicity amplitudes against the squared matrix element derivedfrom an independent computation that directly computes the interference betweenthe tree level and the two loop amplitudes in the HV scheme.• We compare the finite remainder of the squared matrix element obtained from thehelicity amplitudes against the results of [14] at the following physical point sm t = 5 , tm t = − , µ = m t , m t = 1 , (7.3)– 27 – − (cid:15) − (cid:15) − (cid:15) − (cid:15) (cid:98) A (2) , (cid:98) A (2) , − (cid:98) A (2) ,N l ++++ (cid:98) A (2) ,N l +++ − (cid:98) A (2) ,N h ++++ (cid:98) A (2) ,N h +++ − (cid:98) A (2) ,N l ++++ (cid:98) A (2) ,N l +++ − (cid:98) A (2) ,N l N h ++++ (cid:98) A (2) ,N l N h +++ − (cid:98) A (2) ,N h ++++ (cid:98) A (2) ,N h +++ − Table 1 : Mass renormalised two-loop amplitudes (normalised to the tree level amplitude)using the kinematic point in Eq. (7.1).with N l + N h active flavours in the running of α s has been used. In order to extractthe finite remainder with N l + N h active flavours from [14], we subtract the IR polesin the CDR scheme from the renormalised amplitude given in Table 3 of [14]. Asexplained in Section 6, in this paper we have considered evaluation of the iteratedintegrals associated with the curves a and b in the physical region from analyticalresults while leaving the study of the analytic continuation of the integrals whichdepend also on the curve c for the future publication. Therefore, for the evaluationof the amplitude with a single closed top-loop contribution ( A (2) ,N h ) in the physicalregion, we obtain numerical results of the mass renormalised amplitude in the masterintegral representation (Eq. (4.14)). In this representation, the master integrals con-taining elliptic curve c are evaluated using Fiesta [95] and pySecDec [96], while theremaining master integrals are evaluated from their analytical expressions. From thismaster integral representation we further remove the UV and IR poles numerically toobtain the finite remainder.
In this article we have computed a set of analytic helicity amplitudes for the planar cor-rections to top quark pair production via gluon fusion at two loops in QCD. While theseamplitudes have been available for some time in alternative forms, this is the first time thatamplitude level expressions using the massive spinor-helicity formalism have been obtained.We demonstrate that this form is suitable for processing with a rational parametrisationof the kinematics. This has the benefit that the algebra can be performed using finite field– 28 –rithmetic and so combat the growth in algebraic complexity. This method will scale betterthan conventional approaches when considering additional final state particles.Another important new ingredient is the inclusion of top quark loops using analyticexpressions for the master integrals. Firstly we observed that the complete set of masterintegrals contributing to the amplitude were not available in the literature and computedthem using the differential equation method. These integrals give rise to a set of specialfunctions which have recently been presented in terms of iterated integrals over kernelsincluding three elliptic curves. The growth in analytic complexity encountered in this sectoris considerable and we presented a study of the numerical evaluation as well as analyticsimplification leading to complete cancellation of the universal IR and UV poles.We provide a complete set of independent finite remainders in the ancillary files ac-companying the arXiv version of this article. We have also performed a cross check on theingredients in our computation against the available numerical results.Our work leaves a few open questions that motivate further investigation. The numeri-cal evaluation of the iterated integrals is not at the same level of maturity as the commonlyused MPLs. In particular the analytic continuation allowing a stable and efficient evalua-tion over the whole physical scattering region requires further study. We also notice thatthe iterated integrals obey additional functional relations. A complete understanding ofthese identities could and lead to a substantial reduction in the complexity of the finalexpressions and deserves additional study.
Acknowledgments
We thank Matteo Becchetti and Stefan Weinzierl for helpful discussions. This projecthas received funding from the European Union’s Horizon 2020 research and innovationprogrammes
New level of theoretical precision for LHC Run 2 and beyond (grant agree-ment No 683211) and
High precision multi-jet dynamics at the LHC (grant agreementNo 772009). HBH has been partially supported by STFC consolidated HEP theory grantST/T000694/1.
A Renormalisation constants
In this Appendix, we provide renormalisation constants required to derive the renormalisedamplitudes defined in Eqs. (2.7) and (2.8). Their expansion in the bare strong couplingconstant α s are – 29 – m = 1 + (cid:18) α s π (cid:19) δZ (1) m + (cid:18) α s π (cid:19) δZ (2) m + O (cid:0) ( α s ) (cid:1) , (A.1) Z t = 1 + (cid:18) α s π (cid:19) δZ (1) t + (cid:18) α s π (cid:19) δZ (2) t + O (cid:0) ( α s ) (cid:1) , (A.2) Z g = 1 + (cid:18) α s π (cid:19) δZ (1) g + (cid:18) α s π (cid:19) δZ (2) g + O (cid:0) ( α s ) (cid:1) , (A.3) α s = α s (cid:26) (cid:18) α s π (cid:19) δZ (1) α s + (cid:18) α s π (cid:19) δZ (2) α s (cid:27) + O ( α s ) . (A.4)We recompute the following d s dependent mass counterterms in the integral representationthat contribute at leading colour δZ (1) m = N c d s − d s (cid:15) − (cid:15) ) m t I (cid:0) (cid:1) , (A.5) δZ (2) m = N c (cid:26) m t (cid:18) − (cid:15) + 6 (cid:15) − (cid:15) ) − d s − (cid:15) − (cid:15) ) + d s (cid:19) I (cid:0) (cid:1) + 1 m t (1 − (cid:15) ) (cid:18) − − (cid:15) + 12 (cid:15) − (cid:15) )(1 − (cid:15) ) − d s − (cid:15) + 10 (cid:15) − (cid:15) ) + d s − (cid:15) (cid:19) I (cid:0) (cid:1)(cid:27) + N c N l − d s − (cid:15) + 4 d s (cid:15) − (cid:15) ) m t I (cid:0) (cid:1) + N c N h (cid:26) − − (cid:15) )( − − (cid:15) + d s (cid:15) + (cid:15) + d s (cid:15) )4 (cid:15) (1 + (cid:15) )(1 − (cid:15) ) m t I (cid:0) (cid:1) + − − (cid:15) + d s (cid:15) + 3 (cid:15) + d s (cid:15) (cid:15) (1 + (cid:15) ) m t I (cid:0) (cid:1)(cid:27) , (A.6)where the tadpole and on-shell sunrise integrals are defined by I (cid:0) (cid:1) = (cid:0) m t (cid:1) (cid:15) (cid:90) d d k iπ d/ e − (cid:15)γ E k − m t , (A.7) I (cid:0) (cid:1) = (cid:0) m t (cid:1) (cid:15) (cid:90) d d k iπ d/ e − (cid:15)γ E d d k iπ d/ e − (cid:15)γ E k − m t )( k − m t ) , (A.8) I (cid:0) (cid:1) = (cid:0) m t (cid:1) (cid:15) (cid:90) d d k iπ d/ e − (cid:15)γ E d d k iπ d/ e − (cid:15)γ E k − m t ) k ( k + k + p ) , (A.9) I (cid:0) (cid:1) = (cid:0) m t (cid:1) (cid:15) (cid:90) d d k iπ d/ e − (cid:15)γ E d d k iπ d/ e − (cid:15)γ E k − m t )( k − m t )(( k + k + p ) − m t ) , (A.10)where p = m t in Eqs. (A.9) and (A.10). The normalisation of the integrals appearing inEqs. (A.5) and (A.6) is chosen such that it matches the convention adopted in Eq. (2.2). Wehave checked that the integral representation of mass counterterms reproduce the knownintegrated forms given in Ref. [97]. The wavefunction and strong coupling renormalisation– 30 –onstants at one loop are [27] δZ (1) t = C (cid:15) ( m t ) − (cid:15) C F (cid:18) − (cid:15) − − (cid:15) (cid:19) , (A.11) δZ (1) g = C (cid:15) ( m t ) − (cid:15) T F N h (cid:18) − (cid:15) (cid:19) , (A.12) δZ (1) α s = C (cid:15) e − γ E Γ(1 + (cid:15) ) (cid:18) − β (cid:15) (cid:19) , (A.13)and two loops [27, 98, 99] δZ (2) t = C (cid:15) ( m t ) − (cid:15) C F (cid:26) C F (cid:18) (cid:15) + 514 (cid:15) + 4338 − ζ (3) + 96 ζ (2) log(2) − ζ (2) (cid:19) + C A (cid:18) − (cid:15) − (cid:15) − ζ (3) − ζ (2) log(2) + 30 ζ (2) (cid:19) + N l T F (cid:18) (cid:15) + 9 (cid:15) + 8 ζ (2) + 592 (cid:19) + N h T F (cid:18) (cid:15) + 193 (cid:15) + 113918 − ζ (2) (cid:19) + O ( (cid:15) ) (cid:27) , (A.14) δZ (2) g = C (cid:15) ( m t ) − (cid:15) T F N h (cid:26) C F − − (cid:15) + 4 (cid:15) ) (cid:15) (2 − (cid:15) − (cid:15) + 4 (cid:15) )+ C A − − (cid:15) + (cid:15) + 15 (cid:15) − (cid:15) ) (cid:15) (6 + 7 (cid:15) − (cid:15) − (cid:15) + 4 (cid:15) ) + N h T F (cid:15) (cid:27) , (A.15) δZ (2) α s = C (cid:15) e − γ E Γ (1 + (cid:15) ) (cid:18) β (cid:15) − β (cid:15) (cid:19) , (A.16)where C (cid:15) = (4 π ) (cid:15) Γ(1 + (cid:15) ) . The beta function coefficients are β = 113 C A − T F ( N l + N h ) , (A.17) β = 343 C A − C A T F ( N l + N h ) − C F T F ( N l + N h ) . (A.18) B Numerical results for one-loop amplitudes
We present in Table 2 numerical results for one-loop mass-renormalised amplitude at thekinematic point given in Eq. (7.1), evaluated through O ( (cid:15) ) for various fermion loop con-tributions defined in Eq. (2.3). These one-loop amplitudes are required in the computationof the pole terms in Eqs. (4.16) and (4.17). References [1] P. Nason, S. Dawson and R. K. Ellis,
The Total Cross-Section for the Production of HeavyQuarks in Hadronic Collisions , Nucl. Phys. B (1988) 607–633.[2] P. Nason, S. Dawson and R. K. Ellis,
The One Particle Inclusive Differential Cross-Sectionfor Heavy Quark Production in Hadronic Collisions , Nucl. Phys. B (1989) 49–92. – 31 – − (cid:15) − (cid:15) (cid:15) (cid:15) (cid:98) A (1) , -2 -1.5744168 -89.02897164 -497.5712587 -1803.749249 (cid:98) A (1) , − -2 -1.5744168 -9.997841784 -47.34990423 -200.9328154 (cid:98) A (1) ,N l ++++ (cid:98) A (1) ,N l +++ − (cid:98) A (1) ,N h ++++ (cid:98) A (1) ,N h +++ − Table 2 : Mass renormalised one-loop amplitudes (normalised to the tree level amplitude)using the kinematic point in Eq. (7.1). [3] P. Bärnreuther, M. Czakon and A. Mitov,
Percent Level Precision Physics at the Tevatron:First Genuine NNLO QCD Corrections to q ¯ q → t ¯ t + X , Phys. Rev. Lett. (2012) 132001,[ ].[4] M. Czakon and A. Mitov,
NNLO corrections to top-pair production at hadron colliders: theall-fermionic scattering channels , JHEP (2012) 054, [ ].[5] M. Czakon and A. Mitov, NNLO corrections to top pair production at hadron colliders: thequark-gluon reaction , JHEP (2013) 080, [ ].[6] M. Czakon, P. Fiedler and A. Mitov, Total Top-Quark Pair-Production Cross Section atHadron Colliders Through O ( α S ) , Phys. Rev. Lett. (2013) 252004, [ ].[7] M. Czakon, D. Heymes and A. Mitov,
High-precision differential predictions for top-quarkpairs at the LHC , Phys. Rev. Lett. (2016) 082003, [ ].[8] M. Czakon,
A novel subtraction scheme for double-real radiation at NNLO , Phys. Lett. B (2010) 259–268, [ ].[9] R. Bonciani, S. Catani, M. Grazzini, H. Sargsyan and A. Torre,
The q T subtraction methodfor top quark production at hadron colliders , Eur. Phys. J. C (2015) 581, [ ].[10] S. Catani, S. Devoto, M. Grazzini, S. Kallweit, J. Mazzitelli and H. Sargsyan, Top-quark pairhadroproduction at next-to-next-to-leading order in QCD , Phys. Rev. D (2019) 051501,[ ].[11] S. Catani, S. Devoto, M. Grazzini, S. Kallweit and J. Mazzitelli, Top-quark pair production atthe LHC: Fully differential QCD predictions at NNLO , JHEP (2019) 100, [ ].[12] S. Catani, S. Devoto, M. Grazzini, S. Kallweit and J. Mazzitelli, Top-quark pairhadroproduction at NNLO: differential predictions with the
M S mass , JHEP (2020) 027,[ ].[13] M. Czakon, Tops from Light Quarks: Full Mass Dependence at Two-Loops in QCD , Phys.Lett. B (2008) 307–314, [ ].[14] P. Bärnreuther, M. Czakon and P. Fiedler,
Virtual amplitudes and threshold behaviour ofhadronic top-quark pair-production cross sections , JHEP (2014) 078, [ ].[15] L. Chen, M. Czakon and R. Poncelet, Polarized double-virtual amplitudes for heavy-quarkpair production , JHEP (2018) 085, [ ]. – 32 –
16] L. Adams, E. Chaubey and S. Weinzierl,
Simplifying differential equations for multiscalefeynman integrals beyond multiple polylogarithms , Physical Review Letters (Apr, 2017) .[17] L. Adams, E. Chaubey and S. Weinzierl,
Planar Double Box Integral for Top PairProduction with a Closed Top Loop to all orders in the Dimensional RegularizationParameter , Phys. Rev. Lett. (2018) 142001, [ ].[18] L. Adams, E. Chaubey and S. Weinzierl,
Analytic results for the planar double box integralrelevant to top-pair production with a closed top loop , JHEP (2018) 206, [ ].[19] L. Adams and S. Weinzierl, Feynman integrals and iterated integrals of modular forms , 2018.[20] S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow and M. Zeng,
Two-loop integrals forplanar five-point one-mass processes , Journal of High Energy Physics (Nov, 2020) .[21] L. Adams and S. Weinzierl,
The ε -form of the differential equations for feynman integrals inthe elliptic case , Physics Letters B (Jun, 2018) 270?278.[22] C. Bogner, A. Schweitzer and S. Weinzierl,
Analytic continuation and numerical evaluation ofthe kite integral and the equal mass sunrise integral , Nuclear Physics B (2017) 528–550.[23] J. Broedel, C. Duhr, F. Dulat, R. Marzucca, B. Penante and L. Tancredi,
An analyticsolution for the equal-mass banana graph , JHEP (2019) 112, [ ].[24] S. Abreu, M. Becchetti, C. Duhr and R. Marzucca, Three-loop contributions to the ρ parameter and iterated integrals of modular forms , JHEP (2020) 050, [ ].[25] R. Bonciani, A. Ferroglia, T. Gehrmann, D. Maitre and C. Studerus, Two-Loop FermionicCorrections to Heavy-Quark Pair Production: The Quark-Antiquark Channel , JHEP (2008) 129, [ ].[26] R. Bonciani, A. Ferroglia, T. Gehrmann and C. Studerus, Two-Loop Planar Corrections toHeavy-Quark Pair Production in the Quark-Antiquark Channel , JHEP (2009) 067,[ ].[27] R. Bonciani, A. Ferroglia, T. Gehrmann, A. von Manteuffel and C. Studerus, Two-LoopLeading Color Corrections to Heavy-Quark Pair Production in the Gluon Fusion Channel , JHEP (2011) 102, [ ].[28] R. Bonciani, A. Ferroglia, T. Gehrmann, A. von Manteuffel and C. Studerus, Light-quarktwo-loop corrections to heavy-quark pair production in the gluon fusion channel , JHEP (2013) 038, [ ].[29] A. von Manteuffel and C. Studerus, Massive planar and non-planar double box integrals forlight Nf contributions to gg- > tt , JHEP (2013) 037, [ ].[30] S. Di Vita, S. Laporta, P. Mastrolia, A. Primo and U. Schubert, Master integrals for theNNLO virtual corrections to µe scattering in QED: the non-planar graphs , JHEP (2018)016, [ ].[31] M. Becchetti, R. Bonciani, V. Casconi, A. Ferroglia, S. Lavacca and A. von Manteuffel, Master Integrals for the two-loop, non-planar QCD corrections to top-quark pair productionin the quark-annihilation channel , JHEP (2019) 071, [ ].[32] S. Di Vita, T. Gehrmann, S. Laporta, P. Mastrolia, A. Primo and U. Schubert, Masterintegrals for the NNLO virtual corrections to qq → tt scattering in QCD: the non-planargraphs , JHEP (2019) 117, [ ]. – 33 –
33] P. S. Wang,
A p-adic algorithm for univariate partial fractions , in
Proceedings of the FourthACM Symposium on Symbolic and Algebraic Computation , SYMSAC ’81, (New York, NY,USA), pp. 212–217, ACM, 1981. DOI.[34] P. S. Wang, M. J. T. Guy and J. H. Davenport,
P-adic reconstruction of rational numbers , SIGSAM Bull. (May, 1982) 2–3.[35] ISSAC ’06: Proceedings of the 2006 International Symposium on Symbolic and AlgebraicComputation , (New York, NY, USA), ACM, 2006.[36] A. von Manteuffel and R. M. Schabinger,
A novel approach to integration by parts reduction , Phys. Lett.
B744 (2015) 101–104, [ ].[37] T. Peraro,
Scattering amplitudes over finite fields and multivariate functional reconstruction , JHEP (2016) 030, [ ].[38] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflowgraphs , .[39] R. Ellis, W. T. Giele, Z. Kunszt and K. Melnikov, Masses, fermions and generalized D -dimensional unitarity , Nucl. Phys. B (2009) 270–282, [ ].[40] K. Melnikov and M. Schulze,
NLO QCD corrections to top quark pair production and decayat hadron colliders , JHEP (2009) 049, [ ].[41] S. Badger, R. Sattler and V. Yundin, One-Loop Helicity Amplitudes for t ¯ t Production atHadron Colliders , Phys. Rev. D (2011) 074020, [ ].[42] S. Badger, C. Brønnum-Hansen, F. Buciuni and D. O’Connell, A unitarity compatibleapproach to one-loop amplitudes with massive fermions , JHEP (2017) 141, [ ].[43] W. Bernreuther, A. Brandenburg, Z. G. Si and P. Uwer, Top quark spin correlations athadron colliders: Predictions at next-to-leading order QCD , Phys. Rev. Lett. (2001)242002, [ hep-ph/0107086 ].[44] W. Bernreuther, A. Brandenburg, Z. G. Si and P. Uwer, Top quark pair production and decayat hadron colliders , Nucl. Phys. B (2004) 81–137, [ hep-ph/0403035 ].[45] Z. Bern, L. J. Dixon and D. A. Kosower,
Dimensionally regulated pentagon integrals , Nucl.Phys. B (1994) 751–816, [ hep-ph/9306240 ].[46] A. V. Kotikov,
Differential equations method: New technique for massive Feynman diagramscalculation , Phys. Lett. B (1991) 158–164.[47] E. Remiddi,
Differential equations for Feynman graph amplitudes , Nuovo Cim. A (1997)1435–1452, [ hep-th/9711188 ].[48] T. Gehrmann and E. Remiddi,
Differential equations for two loop four point functions , Nucl.Phys. B (2000) 485–518, [ hep-ph/9912329 ].[49] J. M. Henn,
Multiloop integrals in dimensional regularization made simple , Phys. Rev. Lett. (2013) 251601, [ ].[50] S. Catani,
The Singular behavior of QCD amplitudes at two loop order , Phys. Lett.
B427 (1998) 161–171, [ hep-ph/9802439 ].[51] T. Becher and M. Neubert,
Infrared singularities of scattering amplitudes in perturbativeQCD , Phys. Rev. Lett. (2009) 162001, [ ]. – 34 –
52] T. Becher and M. Neubert,
On the Structure of Infrared Singularities of Gauge-TheoryAmplitudes , JHEP (2009) 081, [ ].[53] S. Catani, S. Dittmaier and Z. Trocsanyi, One loop singular behavior of QCD and SUSYQCD amplitudes with massive partons , Phys. Lett. B (2001) 149–160, [ hep-ph/0011222 ].[54] A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang,
Two-loop divergences of scatteringamplitudes with massive partons , Phys. Rev. Lett. (2009) 201601, [ ].[55] A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang,
Two-loop divergences of massivescattering amplitudes in non-abelian gauge theories , JHEP (2009) 062, [ ].[56] R. Kleiss and W. J. Stirling, Spinor Techniques for Calculating p anti-p — > W+- / Z0 +Jets , Nucl. Phys. B (1985) 235–262.[57] A. Hodges,
Eliminating spurious poles from gauge-theoretic amplitudes , JHEP (2013) 135,[ ].[58] G. Rodrigo, Multigluonic scattering amplitudes of heavy quarks , JHEP (2005) 079,[ hep-ph/0508138 ].[59] C. Schwinn and S. Weinzierl, On-shell recursion relations for all Born QCD amplitudes , JHEP (2007) 072, [ hep-ph/0703021 ].[60] N. Arkani-Hamed, T.-C. Huang and Y.-t. Huang, Scattering Amplitudes For All Masses andSpins , .[61] S. Badger, H. Frellesvig and Y. Zhang, A Two-Loop Five-Gluon Helicity Amplitude in QCD , JHEP (2013) 045, [ ].[62] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, M. Kraus, B. Page et al., Caravel: A C++Framework for the Computation of Multi-Loop Amplitudes with Numerical Unitarity , .[63] L. Chen, A prescription for projectors to compute helicity amplitudes in D dimensions , .[64] T. Peraro and L. Tancredi, Physical projectors for multi-leg helicity amplitudes , JHEP (2019) 114, [ ].[65] T. Peraro and L. Tancredi, Tensor decomposition for bosonic and fermionic scatteringamplitudes , .[66] P. Nogueira, Automatic Feynman graph generation , J. Comput. Phys. (1993) 279–289.[67] J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga,
FORM version 4.0 , Comput.Phys. Commun. (2013) 1453–1467, [ ].[68] B. Ruijl, T. Ueda and J. Vermaseren,
FORM version 4.2 , .[69] G. Cullen, M. Koch-Janusz and T. Reiter, Spinney: A Form Library for Helicity Spinors , Comput. Phys. Commun. (2011) 2368–2387, [ ].[70] H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro,
A numerical evaluation ofplanar two-loop helicity amplitudes for a W-boson plus four partons , JHEP (2019) 119,[ ].[71] P. Mastrolia, T. Peraro and A. Primo, Adaptive Integrand Decomposition in parallel andorthogonal space , JHEP (2016) 164, [ ]. – 35 –
72] S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro,
First look at two-loopfive-gluon scattering in QCD , Phys. Rev. Lett. (2018) 092001, [ ].[73] S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro,
Analytic helicity amplitudesfor two-loop five-gluon scattering: the single-minus case , JHEP (2019) 186, [ ].[74] S. Badger, H. B. Hartanto and S. Zoia, Two-loop QCD corrections to
W b ¯ b production athadron colliders , .[75] R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction , .[76] S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations , Int. J. Mod. Phys.
A15 (2000) 5087–5159, [ hep-ph/0102033 ].[77] P. Mastrolia, M. Passera, A. Primo and U. Schubert,
Master integrals for the NNLO virtualcorrections to µe scattering in QED: the planar graphs , JHEP (2017) 198, [ ].[78] L.-B. Chen and J. Wang, Master integrals of a planar double-box family for top-quark pairproduction , Phys. Lett. B (2019) 50–55, [ ].[79] K.-T. Chen,
Iterated path integrals , Bull. Am. Math. Soc. (1977) 831–879.[80] F. Brown, Iterated integrals in quantum field theory , Geometric and Topological Methods forQuantum Field Theory (01, 2011) .[81] A. B. Goncharov,
Multiple polylogarithms, cyclotomy and modular complexes , 2011.[82] S. Moch, P. Uwer and S. Weinzierl,
Nested sums, expansion of transcendental functions, andmultiscale multiloop integrals , Journal of Mathematical Physics (Jun, 2002) 3363?3386.[83] J. M. Borwein, D. M. Bradley, D. J. Broadhurst and P. Lisonek, Special values of multiplepolylogarithms , Trans. Am. Math. Soc. (2001) 907–941, [ math/9910045 ].[84] C. Duhr and F. Dulat,
PolyLogTools — polylogs for the masses , JHEP (2019) 135,[ ].[85] E. Panzer, Feynman integrals and hyperlogarithms . PhD thesis, Humboldt U., 2015. . 10.18452/17157.[86] C. W. Bauer, A. Frink and R. Kreckel,
Introduction to the GiNaC framework for symboliccomputation within the C++ programming language , J. Symb. Comput. (2002) 1–12,[ cs/0004015 ].[87] J. H. Silverman, The Arithmetic of Elliptic Curves . Springer, 2nd ed., 1986.[88] J. Broedel, C. Duhr, F. Dulat and L. Tancredi,
Elliptic polylogarithms and iterated integralson elliptic curves. Part I: general formalism , JHEP (2018) 093, [ ].[89] H. Frellesvig and C. G. Papadopoulos, Cuts of feynman integrals in baikov representation , Journal of High Energy Physics (Apr, 2017) .[90] F. Brown,
Multiple modular values and the relative completion of the fundamental group of m , , 2017.[91] J. Broedel, C. Duhr, F. Dulat and L. Tancredi, Elliptic polylogarithms and iterated integralson elliptic curves II: an application to the sunrise integral , Phys. Rev. D (2018) 116009,[ ].[92] J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, Elliptic symbol calculus: fromelliptic polylogarithms to iterated integrals of Eisenstein series , JHEP (2018) 014,[ ]. – 36 –
93] C. Duhr and L. Tancredi,
Algorithms and tools for iterated Eisenstein integrals , JHEP (2020) 105, [ ].[94] J. Vollinga and S. Weinzierl, Numerical evaluation of multiple polylogarithms , ComputerPhysics Communications (May, 2005) 177?194.[95] A. V. Smirnov,
FIESTA4: Optimized Feynman integral calculations with GPU support , Comput. Phys. Commun. (2016) 189–199, [ ].[96] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk et al., pySecDec: atoolbox for the numerical evaluation of multi-scale integrals , .[97] K. Melnikov and T. van Ritbergen, The Three loop on-shell renormalization of QCD andQED , Nucl. Phys. B (2000) 515–546, [ hep-ph/0005131 ].[98] M. Czakon, A. Mitov and S. Moch,
Heavy-quark production in massless quark scattering attwo loops in QCD , Phys. Lett. B (2007) 147–159, [ ].[99] M. Czakon, A. Mitov and S. Moch,
Heavy-quark production in gluon fusion at two loops inQCD , Nucl. Phys. B (2008) 210–250, [ ].].