Two-loop QCD corrections to Wb\bar{b} production at hadron colliders
CCAVENDISH-HEP-21/01
Two-loop QCD corrections to
W b ¯ b production at hadron colliders Simon Badger, ∗ Heribertus Bayu Hartanto, † and Simone Zoia ‡ Dipartimento di Fisica and Arnold-Regge Center, Universit`a di Torino,and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom (Dated: February 25, 2021)We present an analytic computation of the two-loop QCD corrections to u ¯ d → W + b ¯ b for anon-shell W -boson using the leading colour and massless bottom quark approximations. We performan integration-by-parts reduction of the unpolarised squared matrix element using finite field recon-struction techniques and identify an independent basis of special functions that allows an analyticsubtraction of the infrared and ultraviolet poles. This basis is valid for all planar topologies forfive-particle scattering with an off-shell leg. INTRODUCTION
The production of a W -boson in association with a pairof b -quarks at hadron colliders is of fundamental impor-tance as a background to associated Higgs production.The process is one of a prioritised list of 2 → pp → W + 2 j production and the work presented hererepresents a significant step towards achieving a completeclassification of the missing two-loop amplitudes.The process has been studied extensively at next-to-leading order (NLO) [1–5] and was the first in a set ofoff-shell five-particle amplitudes to be studied using theunitarity method [6, 7]. The present state of the art inphenomenological studies allows full mass effects, showermatching, electro-weak corrections and the inclusion ad-ditional QCD jets [8–10].A numerical computation of the two-loop helicity am-plitudes [11] demonstrated the importance of an efficientanalytic form with a well understood basis of special func-tions. Major steps forward came via efficient numeri-cal evaluation of the differential equations [12] and ana-lytic evaluation in terms the Goncharov Polylogarithms(GPLs) [13, 14]. These results opened the door for a fullyanalytic amplitude computation yet significant challengesremain. The complexity of the external kinematics rep-resents a challenge for integral reduction techniques andthe identification of a minimal basis of special functions isrequired to find analytic simplifications after subtractinguniversal infrared and ultraviolet divergences.Efficient amplitude and integration-by-parts reduction(IBP) [15, 16] using finite field arithmetic [17–27] hasgained significant interest in recent years. Through mul-tiple evaluations of a numerical algorithm [28–31], fullyanalytic forms for planar massless five-particle ampli-tudes have been extracted using a rational parametri-sation of the kinematics [32]. Following a completeunderstanding of a pentagon function basis [33, 34], alarge number of two-loop amplitudes are now available in compact analytic form [35–47]. We have also seenthe first phenomenological predictions at NNLO in QCDfor the production of three photons in hadron collidersafter combination with real-virtual and double real radi-ation [48, 49].In this short letter we outline the extension of thismethod to processes with an additional mass scale. LEADING COLOUR u ¯ d → W + b ¯ b AMPLITUDES
The leading order process consists of two simple Feyn-man diagrams as shown in Fig. 1. We label our processas follows,¯ d ( p ) + u ( p ) → b ( p ) + ¯ b ( p ) + W + ( p ) , (1)where p = p = p = p = 0 and p = m W . The colourdecomposition at leading colour is A ( L ) (1 ¯ d , u , b , ¯ b , W ) = n L g s g W δ ¯ i i δ ¯ i i A ( L ) (1 ¯ d , u , b , ¯ b , W ) , (2)where n = m (cid:15) N c α s / (4 π ) , α s = g s / (4 π ) and m (cid:15) = i (4 π ) (cid:15) e − (cid:15)γ E . g s and g W are the strong and weak cou-pling constants respectively.We interfere the L -loop partial amplitudes A ( L ) inEq. (2) with the tree-level partial amplitude A (0) to ob-tain the unrenormalised L -loop unpolarised squared par-tial amplitude, M ( L ) = (cid:88) spin A (0) ∗ A ( L ) . (3) p p p p p p p p p p FIG. 1. Leading order Feynman diagrams contributing to u ¯ d → W + b ¯ b . a r X i v : . [ h e p - ph ] F e b After the interference with the tree-level amplitude theanalytic expression can be written in terms of scalar in-variants, s = ( p + p ) , s = ( p − p ) , s = ( p + p ) ,s = ( p + p ) , s = ( p − p ) , s = p . (4)Our results are presented after subtraction of infraredand ultraviolet divergences, F ( L ) = M ( L ) − P ( L ) , where P ( L ) takes the well known form [50–53]. The explicitform for our process using the same conventions can befound in Ref. [11]. AMPLITUDE REDUCTION
Feynman diagrams for the u ¯ d → W + b ¯ b scattering aregenerated using Qgraf [54]. In the leading colour ap-proximation, there are 2, 16 and 210 diagrams contribut-ing to the tree level, 1-loop and two-loop amplitudes,respectively. Upon interference of the L -loop partial am-plitude A ( L ) with the tree level partial amplitude A (0) according to Eq. (3), the squared partial amplitude canbe written as M (2) ( { p } ) = (cid:90) (cid:89) i =1 d d k i iπ d/ e − (cid:15)γ E (cid:88) T N T ( d, d s , { k } , { p } ) (cid:81) α ∈ T D α ( { k } , { p } ) , (5)where p i are the external momenta which live in fourdimensions, and k i are the loop momenta in d = 4 − (cid:15) dimensions. We work in the conventional dimensionalregularisation (CDR) scheme where d = d s = g µµ .We treat the γ to be anti-commuting with all other γ matrices in d dimensions and the polarisation sum isperformed in the unitary gauge, i (cid:88) λ ε µ ∗ W ( p , λ ) ε νW ( p , λ ) = − g µν + p µ p ν m W . (6)In obtaining the numerator function N T , the fermiontraces can safely be evaluated in d dimensions since the γ terms from the L -loop numerator cancelled out withthe γ terms from the tree level partial amplitude.To perform the reduction of the amplitude onto a ba-sis of master integrals we first map each topology T toa set of 15 maximal cut or master topologies as shown W + ub ¯ b ¯ d u b ¯ b ¯ dW + ¯ d W + ub ¯ b FIG. 2. Sample Feynman diagrams in the leading colour two-loop u ¯ d → W + b ¯ b amplitude. T
43 1 52 T
32 4 15 T T T T T T T T T
51 2 34 T
14 5 23 T
43 1 52 T
32 4 15 T
25 3 41
FIG. 3. Topologies with maximum number of propagators. in Fig. 3. The master topologies are then defined with aspanning set of 11 propagators and, after tracking shiftsin the loop momentum, the change of variables for eachtopology T can be computed. The resulting squared par-tial amplitude is now written as a linear combination ofscalar integrals I i M (2) ( { p } ) = (cid:88) i c i ( (cid:15), { p } ) I i ( (cid:15), { p } ) . (7)An analytic form of the unreduced squared matrix ele-ment above is derived using a collection of Form [55, 56]and
Mathematica routines. The integrals appearing inEq. (7) are not all independent. Relations between inte-grals I can be found using IBP identities and the squaredamplitude can be written in terms of an independent setof master integrals as follows M (2) ( { p } ) = (cid:88) i d i ( (cid:15), { p } ) MI i ( (cid:15), { p } ) . (8)The reduction to master integral basis is then per-formed within the FiniteFlow framework [23]. We use
LiteRed [57] to generate the IBP relations in
Math-ematica , together with the Laporta algorithm [58] tosolve them numerically over finite fields. We note thatonly master topologies T − T are included in the IBPsystem since the integrals belonging to master topologies T − T can be mapped onto master topologies T − T .The procedure for performing the reduction onto masterintegrals using IBP relations is of course extremely wellknown, the challenge in this example is one of enormousalgebraic complexity. By encoding the problem withina numeric sampling modular arithmetic we are able toefficiently solve the Laporta system with tensor integralranks of up to five, avoiding all large intermediate expres-sions. For planar topologies such as the ones appearinghere the application of syzygy relations [59–61] to opti-mise the IBP reduction would likely lead to a substantialspeed-up in computation time, although in our case itwas not found to be necessary. We did not perform ananalytic reconstruction after completing the set up of thereduction in FiniteFlow graphs. Instead we continuedto map the amplitude onto a basis of special functions.We did, however, perform a gauge invariance check in themaster integral basis: we modified the numerator func-tions by replacing the loop-amplitude polarisation vectorwith p and the tree one with p , and observed that M (2) vanishes. A BASIS OF SPECIAL FUNCTIONS FOR THEFINITE REMAINDER
There are 202 master integrals contributing to the am-plitude, 196 of them are covered by the 3 independentpentabox master integral topologies, while 6 are of one-loop squared type that involve one-loop massive on-shellbubble integral. We choose the canonical bases of masterintegrals constructed in Ref. [12]. They satisfy differen-tial equations (DEs) [62–65] in the canonical form [66], d −→ MI = (cid:15) (cid:88) i =1 a i d log w i −→ MI , (9)where −→ MI is the set of canonical master integrals for anyof the involved topologies, the a i are constant rationalmatrices, while { w i } i =1 is a set of algebraic functions ofthe external kinematics called letters (see Ref. [12] fortheir definition). The alphabet, i.e. the set of all let-ters, is the same for all planar one-mass five-particle in-tegrals up to two loops, whereas the constant matrices a i depend on the topology. In Ref. [12], the authors alsodiscuss a strategy to evaluate the master integrals nu-merically, based on the solution of the DEs (9) in termsof generalised power series [67]. More recently, analyticexpressions of the canonical master integrals in terms ofGPLs [68–70] have become available [13, 14]. Both ap-proaches allow for the numerical evaluation of the mas-ter integrals in any kinematic region and with arbitraryprecision. Both approaches, however, also share certaindrawbacks. Whether we reconstruct the prefactors ofthe (cid:15) -components of the master integrals in Eq. (8) or wemap the latter onto monomials of GPLs, we cannot sub-tract the infrared and ultraviolet poles analytically andreconstruct directly the finite remainder.We overcome these issues by constructing a basis out ofthe (cid:15) -components of the canonical master integrals up toorder (cid:15) . The crucial tool we employ in this constructionare Chen’s iterated integrals [71]. We can define themiteratively through d [ w i , . . . , w i n ] s ( s ) = d log w i n [ w i , . . . , w i n − ] s ( s ) , [ w i , . . . , w i n ] s ( s ) = 0 , (10)where s denotes cumulatively the kinematic invariants, s is an arbitrary boundary point, and the iteration startsfrom [] s ( s ) = 1. The depth n of the iterated integral is called transcendental weight. We refer to the notes [72]for a thorough discussion. All GPLs can be rewrittenin terms of iterated integrals. The latter however of-fer two useful advantages. The first is that – conjec-turally – they implement automatically all the functionalrelations. Once a GPL expression is rewritten in termsof iterated integrals in a given alphabet { w i } , finding thefunctional relations becomes a linear algebra problem,as ‘words’ [ w i , . . . , w i n ] with different letters are linearlyindependent. The second is that it is completely straight-forward to write out the solution of the canonical DEs (9)in terms of iterated integrals. Eq. (9) in fact implies thefollowing differential relation between consecutive com-ponents of the (cid:15) expansion of the master integrals, d −→ MI ( k ) = (cid:88) i =1 a i d log w i −→ MI ( k − , ∀ k ≥ , (11)where −→ MI ( k ) is the O ( (cid:15) k ) term of the master integrals.Comparing Eq. (11) to Eq. (10), we see that the iter-ated integral expressions of −→ MI ( k ) is obtained by addinga letter to the right of those of the previous order, shuf-fling them as prescribed by the constant matrices a i , andadding the boundary values. The master integrals arenormalised to start from O ( (cid:15) ) and so the O ( (cid:15) k ) compo-nents have transcendental weight k .We used the GPL expressions of Refs. [13, 14] to com-pute the values of the master integrals in an arbitrarypoint s with 1100-digit precision. Using the PSLQ algo-rithm [73], we determined the integer relations among theboundary values, and rewrote them in terms of a basis oftranscendental constants. Next, we used the differentialequations provided by Ref. [12] to express the relevantmaster integrals in terms of iterated integrals. This al-lowed us to determine a minimal set of linearly indepen-dent integral components, order by order in (cid:15) up to (cid:15) .We denote these functions by { f ( w ) i } , where w = 1 , . . . , f ( w ) i corresponds to an (cid:15) component of the master integrals, we can evaluate themnumerically using the methods of Refs. [12–14], with theadditional advantage that they are linearly independent.In order to subtract the poles analytically, we needto be able to write in the same basis also the subtrac-tion term. From the transcendental point of view, thelatter is given by the product of certain logarithms andtranscendental constants coming from the anomalous di-mensions – π and ζ – times the one-loop amplitude.In order to accommodate this in the basis, we add thetranscendental constants as elements, and work out therelations between the functions at each weight and prod-ucts of lower-weight ones using the shuffle algebra of theiterated integrals. As a result, the functions in the basis { f ( w ) i } are indecomposable, i.e. they cannot be rewrittenin terms of lower-weight elements of the basis.Armed with this function basis, we can proceed withthe reconstruction of the two-loop finite remainder. Wemap the master integrals appearing in Eq. (8) onto amonomial basis of the functions { f ( w ) i } , which we denoteby { m i ( f ) } , and perform a Laurent expansion in (cid:15) up to O ( (cid:15) ). We do the same for the subtraction term P (2) .The resulting finite remainder, F (2) ( { p } ) = (cid:88) i e i ( { p } ) m i ( f ) + O ( (cid:15) ) , (12)is indeed free of (cid:15) poles. We set s = 1 to simplifythe reconstruction. The dependence can be recovered aposteriori through dimensional analysis. The coefficients e i ( { p } ) in F (2) are not all independent. We find thelinear relations between them and a set of additional co-efficients we supply as ansatz. We used tree-level expres-sions, coefficients from the one-loop amplitude and fromthe unreduced scalar integrals. Through these linear re-lations we rewrite the complicated coefficients in F (2) in terms known coefficients from the ansatz and simplerones, which finally have to be reconstructed. Moreover,we simplify the reconstruction of the remaining coeffi-cients by partial fractioning them with respect to s .First we determine the denominator factors by comput-ing a univariate slice and matching it against an ansatzmade of letters w i . Using the information about the de-nominator and the polynomial degree in the numerator,we construct an ansatz for the partial-fractioned expres-sions of the coefficients. Then we fit the ansatz with anumerical sampling. See Refs. [74, 75] for recent work onmultivariate partial fractioning. To emphasize the effec-tiveness of our strategy, note that the coefficients of thetwo-loop amplitude written in terms of GPL monomialshave maximal degree 62. The maximal degree drops to 54when we use the basis of special functions { f ( w ) i } in thefinite remainder, and then to 31 in the remaining 4 vari-ables after partial fractioning. The reconstruction finallyrequired 38663 sample points over 2 prime fields, gaininga factor of 7 in the reconstruction time with respect tothe GPL-based approach [76]. The reconstructed ana-lytic expressions are further simplified using the Multi-variateApart package [75].The iterated integrals expression of the f ( w ) i functionsallow us to study the analytic structure of the finite re-mainder in a very convenient way. Interestingly, we ob-serve that certain letters do not appear. As it was al-ready noted in Ref. [12], the last 9 letters do show upin any two-loop amplitude up to order (cid:15) . Out of therelevant 49 letters, 6 ( w i with i ∈ { , , , , , } )appear in the master integrals but cancel out in the two-loop amplitude truncated at O ( (cid:15) ). Finally, the letter w = 4 i(cid:15) µνρσ p µ p ν p ρ p σ is present in the two-loop am-plitude, but cancels out in the finite remainder. Thisletter has already been observed to exhibit the same be-haviour in all the known massless two-loop five-particleamplitudes [36–45, 77], which has spawned interest in thecontext of cluster algebras [78]. M (2) /M (0) (cid:15) =0 (cid:15) − (cid:15) − -2.19718713546 (cid:15) − − . . i(cid:15) − − . − . i(cid:15) . − . i TABLE I. Numerical results for the leading colour two-loop squared partial amplitude, M (2) , normalised to the treelevel squared partial amplitude in 4 dimensions, M (0) (cid:15) =0 , atthe physical point { s = 5 , s = − / , s = 11 / , s =17 / , s = − / , s = 1 / } . As regards the numerical evaluation, we propose astrategy based on the generalised power series solutionof the DEs [67] applied not to the master integrals, butdirectly to the basis of special functions. If we rescaleeach special function in the basis f ( w ) i by a power of (cid:15) corresponding to its weight, the ensuing list of functions (cid:126)v = { (cid:15) w f ( w ) i , } satisfies a system of DEs in the canoni-cal form (9). This follows from the differential propertyof the iterated integrals (10). Differently from the DEsfor the master integrals, the DEs for the special func-tions contain only the minimal amount of informationnecessary to evaluate the finite remainder. For instance,instead of evaluating all the weight-4 functions that mayappear in any one-mass two-loop five-particle amplitude,we can restrict ourselves to evaluating only the 19 linearcombinations that actually appear in F (2) . We there-fore define a new basis of special functions, { g ( w ) i } , whichat weight four includes the aforementioned 19 combina-tions of f (4) i ’s, at weight three contains only the f (3) i ’sappearing in F (2) and in the derivatives of the g (4) i ’s,and so on down to weight zero. The resulting DEs aremuch simpler that those for the master integrals. For in-stance, they are by-construction free of the letters whichdo not appear in the finite-remainder, and their dimen-sion is smaller than the number of master integrals for allthe relevant families. Finally, we evaluate the g ( w ) i ’s bysolving the corresponding DEs using the Mathematica package
DiffExp [79]. We compute the boundary val-ues in an arbitrary point in the physical scattering regionthrough the correspondence between the g ( w ) i ’s and themaster integral components.The complete analytic expression of the two-loop fi-nite remainder in terms of rational coefficients and spe-cial functions is included in the ancillary files, togetherwith the differential equation and the boundary valuesnecessary to evaluate the latter numerically. For the con-venience of future cross-checks, we provide the numericalvalues of M (2) at one phase space point in Table I. α s = 0 . F i n i t e R e m a i nd e r [ G e V − ] M (0) + (cid:0) α s π (cid:1) { F (1) } + (cid:0) α s π (cid:1) { F (2) } M (0) + (cid:0) α s π (cid:1) { F (1) } M (0) . . . . . . . . . x R a t i o t o M ( ) FIG. 4. The finite remainder for u ¯ d → W + b ¯ b at one and twoloops as a function of the variable x defined in Eq. (13). DISCUSSION AND OUTLOOK
The results we have obtained represent a major stepforward and open the door to phenomenological applica-tions. The identification of a basis of special functions hasresulted in a substantial speed up over previous studiesas well as uncovering explicit cancellations and reduc-tion in complexity. To demonstrate the suitability forphenomenological applications we present the evaluationon a univariate slice of the physical phase space. For thiswe use a parametrisation in terms of energy fractions andangles of the final state, p = x √ s (1 , , , ,p = x √ s (1 , cos θ, − sin φ sin θ, − cos φ sin θ ) ,p = √ s (1 , , , − p − p , (13)where p and p are taken back-to-back along the z -axiswith a total centre-of-mass energy of s . We have cho-sen p to be produced at an elevation of π from the z -axis and the on-shell phase space conditions imposecos θ = 1 + x x (cid:16) − x − x − m W s (cid:17) . In Fig. 4 we plotvalues of the one- and two-loop finite remainders against x for a configuration with φ = 0 . , m W = 0 . , s = 1and x = 0 .
6. The special functions were evaluated with
DiffExp [79] using rationalized values of the invariants.An average evaluation time of 260s over 1000 points wasobserved and the function is smooth and stable over thewhole region. This demonstrates that even with a basicsetup in
Mathematica a reasonable evaluation time canbe achieved and that realistic phenomenology can now beperformed.The results obtained here pave the way for a broader class of 2 → W -boson andapply equally to the planar sectors of pp → W/Z +2 j (in-cluding decays) and pp → H + 2 j . Going beyond leadingcolour for pp → W/Z + 2 j or any complete pp → H + 2 j amplitudes at two-loops still requires missing informa-tion on the non-planar master integrals, nevertheless webelieve they can be easily incorporated into the strategywe introduce here.We thank Herschel Chawdhry, Thomas Gehrmann,Johannes Henn, Alexander Mitov, Tiziano Peraro andRene Poncelet for numerous insightful discussions. Wealso thank Nikolaos Syrrakos for kindly providing the re-sults of Ref. [14] prior to its publication. This projecthas received funding from the European Union’s Horizon2020 research and innovation programmes New level oftheoretical precision for LHC Run 2 and beyond (grantagreement No 683211),
High precision multi-jet dynam-ics at the LHC (grant agreement No 772009), and
Novel structures in scattering amplitudes (grant agree-ment No 725110). HBH has been partially supportedby STFC consolidated HEP theory grant ST/T000694/1.SZ gratefully acknowledges the computing resources pro-vided by the Max Planck Institute for Physics. ∗ [email protected] † [email protected] ‡ [email protected][1] R. K. Ellis and S. Veseli, Strong radiative corrections toW b anti-b production in p anti-p collisions, Phys. Rev.D , 011501 (1999), arXiv:hep-ph/9810489.[2] F. Febres Cordero, L. Reina, and D. Wackeroth, W- andZ-boson production with a massive bottom-quark pairat the Large Hadron Collider, Phys. Rev. D , 034015(2009), arXiv:0906.1923 [hep-ph].[3] S. Badger, J. M. Campbell, and R. K. Ellis, QCD Correc-tions to the Hadronic Production of a Heavy Quark Pairand a W-Boson Including Decay Correlations, JHEP ,027, arXiv:1011.6647 [hep-ph].[4] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pit-tau, and P. Torrielli, W and Z/γ ∗ boson production inassociation with a bottom-antibottom pair, JHEP ,061, arXiv:1106.6019 [hep-ph].[5] C. Oleari and L. Reina, W +- b ¯ b production inPOWHEG, JHEP , 061, [Erratum: JHEP 11, 040(2011)], arXiv:1105.4488 [hep-ph].[6] Z. Bern, L. J. Dixon, D. A. Kosower, and S. Weinzierl,One loop amplitudes for e+ e- — > anti-q q anti-Q Q,Nucl. Phys. B , 3 (1997), arXiv:hep-ph/9610370.[7] Z. Bern, L. J. Dixon, and D. A. Kosower, One loop am-plitudes for e+ e- to four partons, Nucl. Phys. B , 3(1998), arXiv:hep-ph/9708239.[8] G. Luisoni, C. Oleari, and F. Tramontano, W bbj produc-tion at NLO with POWHEG+MiNLO, JHEP , 161,arXiv:1502.01213 [hep-ph].[9] S. Kallweit, J. M. Lindert, P. Maierh¨ofer, S. Pozzorini, and M. Sch¨onherr, NLO electroweak automation and pre-cise predictions for W+multijet production at the LHC,JHEP , 012, arXiv:1412.5157 [hep-ph].[10] F. R. Anger, F. Febres Cordero, H. Ita, and V. Sotnikov,NLO QCD predictions for W b ¯ b production in associationwith up to three light jets at the LHC, Phys. Rev. D ,036018 (2018), arXiv:1712.05721 [hep-ph].[11] H. B. Hartanto, S. Badger, C. Brønnum-Hansen, andT. Peraro, A numerical evaluation of planar two-loop he-licity amplitudes for a W-boson plus four partons, JHEP , 119, arXiv:1906.11862 [hep-ph].[12] S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow,and M. Zeng, Two-Loop Integrals for Planar Five-PointOne-Mass Processes, JHEP , 117, arXiv:2005.04195[hep-ph].[13] D. D. Canko, C. G. Papadopoulos, and N. Syrrakos,Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg, (2020),arXiv:2009.13917 [hep-ph].[14] N. Syrrakos, Pentagon integrals to arbitrary order in thedimensional regulator, (2020), arXiv:2012.10635 [hep-ph].[15] F. V. Tkachov, A Theorem on Analytical Calculabilityof Four Loop Renormalization Group Functions, Phys.Lett. , 65 (1981).[16] K. G. Chetyrkin and F. V. Tkachov, Integration by Parts:The Algorithm to Calculate beta Functions in 4 Loops,Nucl. Phys. B192 , 159 (1981).[17] P. S. Wang, A p-adic algorithm for univariate partialfractions, in
Proceedings of the Fourth ACM Symposiumon Symbolic and Algebraic Computation , SYMSAC ’81(ACM, New York, NY, USA, 1981) pp. 212–217.[18] P. S. Wang, M. J. T. Guy, and J. H. Davenport, P-adicreconstruction of rational numbers, SIGSAM Bull. , 2(1982).[19] ISSAC ’06: Proceedings of the 2006 International Sym-posium on Symbolic and Algebraic Computation (ACM,New York, NY, USA, 2006) 505060.[20] M. Kauers, Fast Solvers for Dense Linear Systems, Nu-clear Physics B (Proc. Suppl.) , 245 (2008).[21] A. von Manteuffel and R. M. Schabinger, A novel ap-proach to integration by parts reduction, Phys. Lett.
B744 , 101 (2015), arXiv:1406.4513 [hep-ph].[22] T. Peraro, Scattering amplitudes over finite fields andmultivariate functional reconstruction, JHEP , 030,arXiv:1608.01902 [hep-ph].[23] T. Peraro, FiniteFlow: multivariate functional recon-struction using finite fields and dataflow graphs, JHEP , 031, arXiv:1905.08019 [hep-ph].[24] A. V. Smirnov and F. S. Chuharev, FIRE6: FeynmanIntegral REduction with Modular Arithmetic, Comput.Phys. Commun. , 106877 (2020), arXiv:1901.07808[hep-ph].[25] J. Klappert and F. Lange, Reconstructing rational func-tions with FireFly, Comput. Phys. Commun. , 106951(2020), arXiv:1904.00009 [cs.SC].[26] J. Klappert, F. Lange, P. Maierh¨ofer, and J. Usovitsch,Integral Reduction with Kira 2.0 and Finite Field Meth-ods, (2020), arXiv:2008.06494 [hep-ph].[27] F. Caola, A. von Manteuffel, and L. Tancredi, Di-photonamplitudes in three-loop Quantum Chromodynamics,(2020), arXiv:2011.13946 [hep-ph].[28] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, andT. Peraro, First look at two-loop five-gluon scatter- ing in QCD, Phys. Rev. Lett. , 092001 (2018),arXiv:1712.02229 [hep-ph].[29] S. Abreu, F. Febres Cordero, H. Ita, B. Page, andM. Zeng, Planar Two-Loop Five-Gluon Amplitudes fromNumerical Unitarity, Phys. Rev. D , 116014 (2018),arXiv:1712.03946 [hep-ph].[30] S. Badger, C. Brønnum-Hansen, T. Gehrmann, H. B.Hartanto, J. Henn, N. A. Lo Presti, and T. Peraro, Ap-plications of integrand reduction to two-loop five-pointscattering amplitudes in QCD, PoS LL2018 , 006 (2018),arXiv:1807.09709 [hep-ph].[31] S. Abreu, F. Febres Cordero, H. Ita, B. Page, and V. Sot-nikov, Planar Two-Loop Five-Parton Amplitudes fromNumerical Unitarity, JHEP , 116, arXiv:1809.09067[hep-ph].[32] A. Hodges, Eliminating spurious poles from gauge-theoretic amplitudes, JHEP , 135, arXiv:0905.1473[hep-th].[33] T. Gehrmann, J. Henn, and N. Lo Presti, Pentagon func-tions for massless planar scattering amplitudes, JHEP , 103, arXiv:1807.09812 [hep-ph].[34] D. Chicherin and V. Sotnikov, Pentagon Functions forScattering of Five Massless Particles, JHEP , 167,arXiv:2009.07803 [hep-ph].[35] T. Gehrmann, J. Henn, and N. Lo Presti, Analyticform of the two-loop planar five-gluon all-plus-helicityamplitude in QCD, Phys. Rev. Lett. , 062001(2016), [Erratum: Phys.Rev.Lett. 116, 189903 (2016)],arXiv:1511.05409 [hep-ph].[36] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, andT. Peraro, Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case, JHEP , 186,arXiv:1811.11699 [hep-ph].[37] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, andB. Page, Analytic Form of Planar Two-Loop Five-GluonScattering Amplitudes in QCD, Phys. Rev. Lett. ,082002 (2019), arXiv:1812.04586 [hep-ph].[38] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page,and V. Sotnikov, Analytic Form of the Planar Two-LoopFive-Parton Scattering Amplitudes in QCD, JHEP ,084, arXiv:1904.00945 [hep-ph].[39] S. Abreu, B. Page, E. Pascual, and V. Sotnikov,Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders, (2020),arXiv:2010.15834 [hep-ph].[40] H. A. Chawdhry, M. Czakon, A. Mitov, and R. Poncelet,Two-loop leading-color helicity amplitudes for three-photon production at the LHC, (2020), arXiv:2012.13553[hep-ph].[41] S. Abreu, L. J. Dixon, E. Herrmann, B. Page, andM. Zeng, The two-loop five-point amplitude in N = 4super-Yang-Mills theory, Phys. Rev. Lett. , 121603(2019), arXiv:1812.08941 [hep-th].[42] D. Chicherin, T. Gehrmann, J. Henn, P. Wasser,Y. Zhang, and S. Zoia, Analytic result for a two-loop five-particle amplitude, Phys. Rev. Lett. , 121602 (2019),arXiv:1812.11057 [hep-th].[43] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser,Y. Zhang, and S. Zoia, The two-loop five-particleamplitude in N = 8 supergravity, JHEP , 115,arXiv:1901.05932 [hep-th].[44] S. Abreu, L. J. Dixon, E. Herrmann, B. Page, andM. Zeng, The two-loop five-point amplitude in N = 8supergravity, JHEP , 123, arXiv:1901.08563 [hep-th]. [45] S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich,J. Henn, T. Peraro, P. Wasser, Y. Zhang, and S. Zoia,Analytic form of the full two-loop five-gluon all-plus he-licity amplitude, Phys. Rev. Lett. , 071601 (2019),arXiv:1905.03733 [hep-ph].[46] B. Agarwal, F. Buccioni, A. von Manteuffel, and L. Tan-credi, Two-loop leading colour QCD corrections to q ¯ q → γγg and qg → γγq , (2021), arXiv:2102.01820 [hep-ph].[47] G. De Laurentis and D. Maˆıtre, Two-Loop Five-PartonLeading-Colour Finite Remainders in the Spinor-HelicityFormalism, JHEP , 016, arXiv:2010.14525 [hep-ph].[48] H. A. Chawdhry, M. L. Czakon, A. Mitov, and R. Pon-celet, NNLO QCD corrections to three-photon produc-tion at the LHC, JHEP , 057, arXiv:1911.00479 [hep-ph].[49] S. Kallweit, V. Sotnikov, and M. Wiesemann, Tripho-ton production at hadron colliders in NNLO QCD10.1016/j.physletb.2020.136013 (2020), arXiv:2010.04681[hep-ph].[50] S. Catani, The Singular behavior of QCD amplitudes attwo loop order, Phys. Lett. B427 , 161 (1998), arXiv:hep-ph/9802439 [hep-ph].[51] T. Becher and M. Neubert, On the Structure of In-frared Singularities of Gauge-Theory Amplitudes, JHEP , 081, [Erratum: JHEP11,024(2013)], arXiv:0903.1126[hep-ph].[52] T. Becher and M. Neubert, Infrared singularities ofscattering amplitudes in perturbative QCD, Phys. Rev.Lett. , 162001 (2009), [Erratum: Phys. Rev.Lett.111,no.19,199905(2013)], arXiv:0901.0722 [hep-ph].[53] E. Gardi and L. Magnea, Factorization constraints forsoft anomalous dimensions in QCD scattering ampli-tudes, JHEP , 079, arXiv:0901.1091 [hep-ph].[54] P. Nogueira, Automatic Feynman graph generation, J.Comput. Phys. , 279 (1993).[55] J. Kuipers, T. Ueda, J. A. M. Vermaseren, andJ. Vollinga, FORM version 4.0, Comput. Phys. Commun. , 1453 (2013), arXiv:1203.6543 [cs.SC].[56] B. Ruijl, T. Ueda, and J. Vermaseren, FORM version 4.2,(2017), arXiv:1707.06453 [hep-ph].[57] R. N. Lee, Presenting LiteRed: a tool for the Loop In-TEgrals REDuction, (2012), arXiv:1212.2685 [hep-ph].[58] S. Laporta, High precision calculation of multiloop Feyn-man integrals by difference equations, Int. J. Mod. Phys. A15 , 5087 (2000), arXiv:hep-ph/0102033 [hep-ph].[59] J. Gluza, K. Kajda, and D. A. Kosower, Towards a Basisfor Planar Two-Loop Integrals, Phys. Rev. D , 045012(2011), arXiv:1009.0472 [hep-th].[60] H. Ita, Two-loop Integrand Decomposition into MasterIntegrals and Surface Terms, Phys. Rev. D , 116015(2016), arXiv:1510.05626 [hep-th].[61] K. J. Larsen and Y. Zhang, Integration-by-parts reduc-tions from unitarity cuts and algebraic geometry, Phys. Rev. D , 041701 (2016), arXiv:1511.01071 [hep-th].[62] A. V. Kotikov, Differential equations method: New tech-nique for massive Feynman diagrams calculation, Phys.Lett. B , 158 (1991).[63] Z. Bern, L. J. Dixon, and D. A. Kosower, Dimension-ally regulated pentagon integrals, Nucl. Phys. B , 751(1994), arXiv:hep-ph/9306240.[64] E. Remiddi, Differential equations for Feynman graphamplitudes, Nuovo Cim. A , 1435 (1997), arXiv:hep-th/9711188.[65] T. Gehrmann and E. Remiddi, Differential equations fortwo loop four point functions, Nucl. Phys. B , 485(2000), arXiv:hep-ph/9912329.[66] J. M. Henn, Multiloop integrals in dimensional regu-larization made simple, Phys. Rev. Lett. , 251601(2013), arXiv:1304.1806 [hep-th].[67] F. Moriello, Generalised power series expansions for theelliptic planar families of Higgs + jet production at twoloops, JHEP , 150, arXiv:1907.13234 [hep-ph].[68] A. B. Goncharov, Multiple polylogarithms, cyclotomyand modular complexes, Math. Res. Lett. , 497 (1998),arXiv:1105.2076 [math.AG].[69] E. Remiddi and J. A. M. Vermaseren, Harmonic polylog-arithms, Int. J. Mod. Phys. A , 725 (2000), arXiv:hep-ph/9905237.[70] A. B. Goncharov, Multiple polylogarithms and mixedTate motives, (2001), arXiv:math/0103059.[71] K.-T. Chen, Iterated path integrals, Bull. Am. Math. Soc. , 831 (1977).[72] F. Brown, Iterated integrals in quantum field theory, in (2013) pp. 188–240.[73] H. R. P. Ferguson and D. H. Bailey, A Polynomial Time,Numerically Stable Integer Relation Algorithm, RNRTechnical Report RNR-91-032 (1992).[74] J. Boehm, M. Wittmann, Z. Wu, Y. Xu, and Y. Zhang,IBP reduction coefficients made simple, JHEP , 054,arXiv:2008.13194 [hep-ph].[75] M. Heller and A. von Manteuffel, MultivariateApart:Generalized Partial Fractions, (2021), arXiv:2101.08283[cs.SC].[76] Since estimates of evaluation time rely on system specificparameters, we have taken a conservative value for thespeed improvement.[77] S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang, andS. Zoia, Multi-Regge Limit of the Two-Loop Five-PointAmplitudes in N = 4 Super Yang-Mills and N = 8 Su-pergravity, JHEP10