Analytic Form of the Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD

We present the analytic form of the two-loop five-gluon scattering amplitudes in QCD for a complete set of independent helicity configurations of external gluons. These include the first analytic results for five-point two-loop amplitudes relevant for the computation of next-to-next-to-leading-order QCD corrections at hadron colliders. The results were obtained by reconstructing analytic expressions from numerical evaluations. The complexity of the computation is reduced by exploiting physical and analytical properties of the amplitudes, employing a minimal basis of so-called pentagon functions that have recently been classified.

Scattering amplitudes in QCD are beautiful mathematical objects that condense the physics of particle collisions. They play a central role in providing theoretical predictions for experimental observables, such as those measured in particle collisions at the Large Hadron Collider at CERN, making higher-order quantum corrections highly desirable. The higher the corrections one considers, the more intricate the amplitudes become, needing to accommodate richer unitarity and factorization properties. Nevertheless, the analytic properties are so constraining that final results appear very compact in spite of the large computational effort required to obtain them. In this work, we exploit these properties to directly access the compact analytic results from numerical evaluations.
Analytic expressions for two-loop four-gluon amplitudes have been known for more than a decade [1,2]. At higher multiplicity, analytic results have been obtained for all-plus helicity amplitudes up to seven points [3][4][5][6][7][8] and up to five points for the single-minus helicity configuration [9], reaching a complexity where traditional methods start to break down. Conventional approaches to analytic calculations have seen considerable progress [10,11] but involve large intermediate expressions, the size of which is highly sensitive to details of the approach and obscures the simplicity of the result. To sidestep this issue and make the most of recent progress on five-point two-loop integral calculations [12][13][14][15], two groups have used numerical approaches to compute fiveparton two-loop QCD scattering amplitudes [16][17][18][19]. Numerical approaches are more resilient to the swelling of intermediate expressions, and can be improved with finitefield calculations to provide exact results [20]. While numerical results are sufficient for Monte Carlo integration, analytic results increase the evaluation efficiency and facilitate the study of the mathematical structure of amplitudes. It has been pointed out that the analytic form of rational functions in scattering amplitudes can be reconstructed from numerical samples [20,21] as already demonstrated for the two-loop four-point amplitudes [22] and recently for the single-minus five-gluon amplitude [9].
Combined with the better understanding of the analytic properties of the functions appearing in loop amplitudes that was developed over the last few years, the reconstruction of analytic expressions from numerical samples circumvents the prohibitively large intermediate analytic steps and directly targets the simpler final results.
In this letter, we use these techniques to compute a complete set of analytic two-loop five-gluon amplitudes in the limit of a large number of colors N c . We apply the extension of the numerical variant [23][24][25][26] of the unitarity method [27][28][29] to multi-loop computations [22,30,31] to obtain numerical samples on finite fields based on [17,19]. From these samples we obtain analytic expressions for the coefficients in a decomposition of the amplitude in terms of so-called pentagon functions [13]. The efficiency of the reconstruction is increased by finding convenient bases in the space of pentagon functions and by an a priori knowledge of the coefficients of the denominators.
We present for the first time analytic results for the two-loop five-gluon amplitudes known as MHV amplitudes, which are an important ingredient in the nextto-next-to-leading-order (NNLO) QCD corrections for three-jet production at hadron colliders. The results will serve as a testing ground for exploring the structure of two-loop scattering amplitudes in quantum field theory.

SCATTERING AMPLITUDES
The main results of this letter are analytic expressions for two-loop five-gluon scattering amplitudes, in the leading-color approximation of pure Yang-Mills theory. We write the perturbative expansion of a bare five-gluon helicity amplitude as Tr (T a σ(1) T a σ(2) T a σ(3) T a σ(4) T a σ(5) ) where λ = N c g 2 0 /(4π) 2 and g 0 is the bare QCD coupling. The set S 5 /Z 5 denotes all non-cyclic permutations of five indices. The planar amplitudes A (k) are functions of the momenta p σ(i) and helicities h σ(i) . In the leading-color approximation there is a single color structure and it is sufficient to specify a helicity assignment for the ordered set of legs, which we specify by subscripts when required, assuming an ordered set of momentum assignments {p 1 , p 2 , p 3 , p 4 , p 5 }. In this section we describe our approach to the calculation of the A (k) in dimensional regularization with D = 4 − 2ǫ. Although we focus on the leading-color contributions, the approach we use is fully generic and applicable beyond this approximation.
Our calculation of the amplitudes A (k) is performed in the framework of two-loop numerical unitarity [17,22,30,31]. The starting point is the general parametrization of the integrand of the amplitude [30] A (2) with M Γ a set of master integrands, S Γ a set of surface terms, P Γ the set of propagators ρ j associated with propagator structure Γ, and ℓ l the set of loop momenta. The coefficients c Γ,i are determined by exploiting the factorization properties of A(ℓ l ) in loop-momenta limits ℓ l → ℓ Γ l where propagators go on-shell [41]: .
The set T Γ labels all tree amplitudes corresponding to the vertices in the diagram Γ. The sum runs over the Γ ′ with more propagators. The state sum runs over the (scheme dependent) gluon-helicity states of each internal line of Γ. Solving the linear system in Eqn. (3) directly yields the coefficients of master integrals and no further integral reduction is needed. In numerical unitarity, the linear system (3) is constructed and solved numerically. Using finite-field arithmetic removes any issues related to loss of precision in numerical operations, at the price of minor modifications to the standard numerical unitarity approach. We describe our algorithm in [17].

NUMERICAL RECONSTRUCTION
We aim to compute analytic expressions for five-gluon amplitudes by reconstructing them from numerical samples. For a suitably chosen parametrization of phase space, such as momentum twistors [32], the coefficients c Γ,i of Eqn. (2) are rational functions. For concreteness, we use the parametrization [21] where s ij = (p i + p j ) 2 with the indices defined cyclically. The difficulty of analytic reconstruction is determined by the complexity of the functions under study. More precisely, this means that one should attempt to reconstruct rational functions with low total degree of numerator and denominator polynomials. With this aim in mind, we first exploit a series of known physical and analytical properties of the amplitudes, and only apply the reconstruction algorithms to simpler objects that we cannot further constrain. In this section we summarize our approach. The divergence structure of scattering amplitudes is governed by known universal functions [33][34][35][36]. It is thus sufficient to compute the so-called finite remainder [33], defined by subtracting contributions that are related to tree and one-loop amplitudes from a two-loop amplitude. There is freedom in how the remainders are defined, so we now give our definitions. For helicity amplitudes which vanish at tree-level, A (k) (5) TheĀ (k) denote amplitudes normalized to remove any ambiguity related to overall phases. In the case of amplitudes that vanish at tree level, we normalize to the leading order in ǫ of the (finite) one-loop amplitude. For the MHV amplitudes, A (k) −∓±++ , which we normalize to the corresponding tree amplitude, we define whereβ i are the coefficients of the QCD β-function divided by N c and S ǫ = (4π) ǫ e −ǫγE . I (1) and I (2) are the standard Catani operators at leading color, and Eqn. (6) is obtained by writing Catani's prediction for the poles of a bare amplitude. Precise expressions for the operators in our conventions can be found in Appendix B of Ref. [19]. We note that in both Eqns. (5) and (6) we require one-loop amplitudes expanded up to order ǫ 2 . The amplitudes we are concerned with in this letter can be written as linear combinations of special functions, the so-called multiple polylogarithms (MPLs). The coefficients in this linear combination are functions of the external data and of the dimensional regulator ǫ, just as in the decomposition of Eqn. (2). It is by now well understood that there are relations between different MPLs, and newly developed mathematical tools allow to find such relations in an algorithmic way [37][38][39]. It is then possible to construct a basis for the space of special functions required for planar five-point massless amplitudes up to two-loops, and this was done in [13] where a set of pentagon functions was defined. Amplitudes can then be written asĀ where x = {x 0 , x 1 , x 2 , x 3 , x 4 } and thec k,i ( x) are rational functions of the twistor variables. B denotes the basis of pentagon functions h i . There is a notion of (transcendental) weight associated with MPLs, and thus with pentagon functions, that allows to organize the basis B.
For instance, there is a single element of weight 0, the trivial function 1. At weight one, there are five elements, ln(−s i,i+1 ). At weight two, there are products of weightone functions, irreducible weight-two functions, and a new constant, π 2 . The same pattern holds at weight three and four, the highest weight required for two-loop amplitudes at order ǫ 0 . To make the simplifications that one expects to find in the remainders explicit, we write both one-and two-loop amplitudes in terms of pentagon functions so that the remainders themselves are written as combinations of pentagon functions: The coefficients r i ( x) have lower total degrees compared with thec k,i ( x) of Eqn. (7), but to further increase the efficiency of the reconstruction we find it useful to implement a series of improvements by investigating the structure of the coefficients on generic 'univariate slices'. Such slices are defined by parametrizing the twistor variables x by a single variable t, x i = a i + b i t, with the a i and b i chosen randomly in the finite field (for high enough cardinality, this ensures the x i do not satisfy simple relations leading to artificial simplifications). The coefficients r i ( x) are themselves univariate rational functions of t, r i ( x(t)) = r i (t). Importantly, on such univariate slices the degrees of the numerator and denominator of the rational functions r i (t) correspond to the total degrees of r i ( x) in the x j . We use these slices to probe the complexity of the functions and find simplifications. Firstly, we classify the pole structures of the coefficient functions r i ( x). A similar classification has been exploited for the computation of one-loop QCD amplitudes in the past [26]. On physical grounds we expect the pole structure to be determined by the so-called alphabet of the pentagon functions. The alphabet determines the points in phase space where the pentagon functions (or their discontinuities) have logarithmic singularities, and they provide a natural ansatz for the poles of the coefficients. The five-point planar alphabet can be written in terms of 26 letters W i , see Ref. [13], which we rewrite in terms of twistor variables. This gives a set A of 26 independent letters A = {w i ( x)}. Indeed, we find that the denominators do factorize into products of letters, The integer exponents q ij are determined by matching this ansatz on univariate slices. The computation of the r i ( x) is then reduced to the much simpler multi-variate polynomial reconstruction of the numerators n i ( x). Secondly, it is expected that cancellations between different basis functions in B take place in exceptional kinematic configurations, which implies relations between their coefficients. This motivates a search for a different basis of pentagon functions with coefficients of lower total degree. To find this new (helicity-dependent) basis, we construct linear combinations of coefficients and solve for phase space independent a i,k such that the numerators N k ( x, a i,k ) factorize a subset of the w j ∈ A. This can be performed on univariate slices by only accepting solutions which are invariant over multiple slices. The matrix a i,k allows to change to a new basis B ′ in the space of special functions, in which remainders can be decomposed as in Eqn. (8), with coefficients r ′ i ( x) whose numerators n ′ i ( x) are polynomials of lower total degree than those of Eqn. (9).
We believe that a deeper understanding of these two improvements might be fruitfully exploited in the future to tackle more complex computations. In particular, it would be interesting to better understand how the branch structure of pentagon functions is related to the poles of the coefficients, and if this suggests a superior basis of pentagon functions with simpler coefficients.

IMPLEMENTATION
Let us now discuss how we applied the techniques of the previous section to perform the analytic reconstruction. We first compute numerical values for the coefficients in a decomposition of one-and two-loop amplitudes in terms of master integrals, using our C++ implementation of numerical unitarity on a finite field [17,19]. We then introduce expressions for the master integrals in terms of pentagon functions and compute the remainders in Eqn. (5) or (6). After changing the basis of pentagon helicityc k,i (t) ri(t) n ′ i (t) wj's in denominator + + + + + t 34 /t 28 t 10 /t 4 t 10 3 − + + + + t 50 /t 42 t 35 /t 28 t 35 14 − − + + + t 70 /t 65 t 50 /t 45 t 40 17 − + − + + t 84 /t 82 t 68 /t 66 t 53 20 functions and multiplying by the predetermined denominator structures, we arrive at an evaluation of the numerator polynomials n ′ i ( x). Table I summarizes the impact of the improvements, the difficulty in the reconstruction of each amplitude, and the number of different letters appearing in the denominators of each helicity. Given the simplicity of R ±++++ , we did not implement any basis change for these helicities. We observe a drastic increase in complexity from top to bottom, as expected from oneloop amplitudes [40]. Comparing the second and third columns we find that the coefficients in Eqn. (8) are indeed simpler than those in Eqn. (7). We also show that not all letters contribute to the denominator in Eqn. (9), even for the most complicated helicity.
The reconstruction of the multivariate polynomials is done with an in-house implementation of the algorithm presented in Ref. [21]. For a detailed description we refer to the original article. The algorithm [21] correlates the sampled phase space to the function it reconstructs, learning about simplifications as it progresses. For this reason, it is not trivial to parallelize. In our implementation this is addressed by anticipating a suitable superset of the points for a set of functions. This implementation is well adapted to evaluations on modern computer clusters.
In order to demonstrate the efficiency of our approach we now summarize the computational requirements for the numerical reconstruction of the presented amplitudes. The dense nature of the algorithm implies that the complexity grows as the number of terms in a polynomial of total degree R in n variables, R+n n . The functions we reconstruct are dimensionless and thus only depend on the four variables x 0 , x 1 , x 2 and x 3 , see Eqn. (4). For example, for a generic polynomial with R = 53, we would need around 400,000 evaluations. In practice, we observe that the polynomials we reconstruct are not completely generic and we require fewer evaluations. For the most complicated remainder, R −+−++ , we require 237,098 phase-space points for the reconstruction of all the functions n ′ i (x). For R −−+++ , we require 85,979 phase-space points. The average time for the numerical computation of the remainder of an MHV amplitude in a finite-field is 4 minutes per phase-space point. The reconstruction of the remainders R ±++++ is simpler. We perform the computation on a finite field of cardinality O(2 31 ). For the MHV amplitudes we also evaluate on a second finite field and apply the Chinese remainder theorem. With this, we are able to rationally reconstruct the result.

RESULTS
We have validated the approach described above by computing the one-loop amplitudes to all orders in ǫ (extending the results of [40]). We then used it to obtain the analytic form of all planar five-gluon two-loop amplitudes in the leading-color approximation, in the 't Hooft-Veltman regularization scheme. Our results are included in a set of ancillary files [42] which we now describe. The files contain analytic expressions for the remainders of four helicity configurations-A ±++++ , see Eqn. (5), and A −∓±++ , see Eqn. (6)-from which any other helicity configuration can be obtained by symmetry or permutation of the momenta. These remainders are written as a linear combination of pentagon functions, similar to Eqn. (8) but with the adjusted bases. This allows to use the library provided in [13] to evaluate the remainders. We also include analytic expressions for the one-loop amplitudes, which are required to obtain the two-loop amplitudes from their remainder. These are written as a linear combination of master integrals in the one-loop basis of Ref. [25] and valid to all orders in ǫ. We provide expressions for the expansion of the one-loop master integrals through order ǫ 2 , written in terms of pentagon functions, so that they can easily be combined with the expressions for the remainders. Finally, we include a script that assembles all components to evaluate a two-loop amplitude.
The results we present are valid in the Euclidean region. Since they are written in terms of pentagon functions [13] there is minimal work to extend the results to all kinematic regions. We checked our expressions against recently computed two-loop amplitudes [3,5,7,9,16,17,19]. In the ancillary files we include numerical values for the pentagon functions at the phasespace point of Ref. [19] and a script that numerically evaluates the amplitudes at that point.

CONCLUSIONS
We have presented the analytic form of a complete set of planar two-loop five-gluon amplitudes in the leadingcolor approximation. The amplitudes were obtained by reconstructing the analytic expressions from numerical samples on finite fields, computed in the framework of two-loop numerical unitarity. This allows to make the most of the resilience of numerical calculations to handle intermediate steps and to only target the analytic form of final expressions, which is constrained by various physical properties. By focusing on finite remainders, we reduced the complexity of the objects to reconstruct. Efficiency of the reconstruction was further enhanced by an a priori determination of the denominators and changes of basis that considerably reduce the total degrees of the numerators of the coefficients. Through this process, we reduced the calculation of the coefficients to the reconstruction of a multi-variate polynomial of relatively low total degree, rendering the most complex MHV amplitudes easily reconstructible on a modern computer cluster.
Our results include the first computation of the analytic form of the five-gluon MHV amplitudes at two loops. They are an important contribution to the NNLO QCD corrections to three-jet production at hadron colliders. While in principle the numerical evaluation of amplitudes is sufficient, the efficiency requirements for phase-space integration over the final states are high due to helicity and color sums. The analytic expressions that we provide will help to control both the evaluation times and the numerical stability of future phenomenological studies.
We expect our computational approach to greatly contribute to formal developments in the study of scattering amplitudes in quantum field theory, and to the new era of precision QCD in high-energy physics.