Double-Virtual NNLO QCD Corrections for Five-Parton Scattering: The Gluon Channel

We compute the two-loop helicity amplitudes for the scattering of five gluons, including all contributions beyond the leading-color approximation. The analytic expressions are represented as linear combinations of transcendental functions with rational coefficients, which we reconstruct from finite-field samples obtained with the numerical unitarity method. Guided by the requirement of removing unphysical singularities, we find a remarkably compact generating set of rational coefficients, which we are able to display entirely in the manuscript. We implement our results in a public code, which provides efficient and reliable numerical evaluations for phenomenological applications.


I. INTRODUCTION
The long-term perspective of the Large Hadron Collider (LHC) at CERN serves as a compelling reason to explore new ways to advance our understanding of high-energy particle collisions beyond the level of detail and precision currently attainable.
Among the central processes under study through hadron collisions is the production of multiple jets.Notably, the recent impressive measurements of the strong coupling constant at high momentum transfer [1][2][3] crucially depend upon the cutting-edge nextto-next-to-leading order (NNLO) QCD predictions for three-jet production [1] (see also [4]).These predictions relied upon the leading-color approximation for double-virtual corrections [5], contributing on average about 10% [1].This highlights the potential importance of including subleading-color effects.Moreover, the observation that the subleading-color effects can be enhanced in certain differential observables [6] underscores the necessity to study three-jet production in full color.
Looking towards the future, five-parton scattering at two loops is also a crucial ingredient in advancing towards the N 3 LO precision frontier for di-jet production in hadron collisions.A related intriguing application lies in the explicit examination of the breakdown of collinear factorization in QCD, which may occur at the third order in perturbation theory when subleading color effects are taken into account [7][8][9].
In this letter we focus on the five-gluon channel and derive compact analytic expressions for all twoloop five-gluon helicity amplitudes, including for the first time all contributions beyond the leading-color approximation.The remaining quark channels will be presented in the followup publication [10].
In this work we build upon this remarkable progress and, in particular, leverage the computational framework established in refs.[31][32][33][34][35].We further delve into a limiting aspect of the analytic reconstruction methods in amplitude computations.Conventional approaches, in essence, employ generic rational ansätze that involve numerous unphysical singularities and redundant parameters.In stark contrast amplitudes often assume concise representations in partial fractioned form.This leads to the intriguing question: can analytic reconstruction directly yield the compact results?While we delay a more thorough discussion of this question to future work [36], we present evidence that a positive outcome can be achieved by leveraging information about the functions' singularities and residues.In fact, considering two-loop five-gluon scattering in full color, we obtain a representation that seamlessly fits into the appendix of this letter.With the exception of special helicity configurations [37][38][39], this level of simplicity has been notably elusive for fivepoint two-loop scattering.

arXiv:2311.10086v3 [hep-ph] 14 May 2024
Finally, we provide a C++ library for fast numerical evaluation of the NNLO hard function that is ready for use in cross-section computations.Together with the upcoming results for the quark channels that we will make available in the followup publication [10], these results will provide crucial input for NNLO cross-section computations.
Note added: while this work was in preparation, we became aware of ref. [40], which reports partially overlapping results.We thank the authors for the numerical comparison of our results and for coordinating the publications.

II. NOTATION AND CONVENTIONS
We consider the O α 2 s corrections to the scattering of five-gluons.This requires the computation of two-loop five-gluon scattering amplitudes, which we obtain omitting the contributions from the massive top quark.Furthermore, we treat all quarks as massless states.The contributing partonic process is 3 ) + g(p h4 4 ) + g(p h5 5 ) .
Here p i and h i denote the momentum and the helicity of the i th particle, respectively.Unless stated otherwise, throughout this paper, momenta and helicity labels are understood in the all-outgoing convention.Representative Feynman diagrams for the two-loop contributions are shown in fig. 1.
Figure 1: Representative Feynman diagrams for two-loop five-gluon amplitudes.Solid lines represent closed massless quark loops.

A. Kinematics
The process involves five massless particles.The underlying scattering kinematic can therefore be specified by five Mandelstam invariants {s 12 , s 23 , s 34 , s 45 , s 15 }, as well as the parity-odd contraction of four momenta tr 5 = tr(γ 5 / p 1 / p 2 / p 3 / p 4 ).To represent the dependence of scattering amplitudes on the particles' helicities we use twocomponent spinors, λ α i and λ α i , with i ∈ {1, . . ., 5}.We define the invariant contractions of spinors as which are related to the Mandelstam invariants through s ij = ⟨ij⟩[ji] (see e.g.[41] for matching conventions).We will also use longer spinor contractions, in particular Finally, we can express tr 5 as a polynomial in spinor brackets as A little-group transformation of the i th leg with helicity h i reads (λ i , λi ) → (z i λ i , λi /z i ).Under this transformation, the helicity amplitudes transform as A → z −2hi i A. We refer to the exponent of the z i as the little-group weight.

B. Color space
The external gluons are in the adjoint representation of SU (N c ) and carry indices ⃗ a = {a 1 , . . ., a 5 } which run over N 2 c − 1 values.We explicitly represent the five-gluon amplitudes in the color space through the trace basis as [42] A where tr(i 1 , . . ., i n ) = tr(T ai 1 • • • T ai n ), and T ai are the hermitian and traceless generators of fundamental representation of SU (N c ).The permutation σ = {i 1 , . . ., i 5 } acts on all external-particle labels as σ(i) = σ i .The sums run over all permutations that do not leave the respective traces invariant.Thus, the first sum runs over 24 elements, while the second one runs over 20 elements.
The generators T a are normalized as, and fulfill the commutator relations, The amplitudes A i admit an expansion in terms of the bare QCD coupling constant α 0 s = (g 0 s )2 /(4π), with L denoting the number of loops.In this work we consider an arbitrary number of light quark flavors in the loops that is denoted by N f .The fivegluon amplitudes can be further expanded in powers of N c and N f through two loops as follows, A (2,0) + A (2),(0,0) (10d) (1,0)  (10e) 2) .
The coefficients A (L),(nc,n f ) , which we call partial amplitudes, are uniquely identified by the three integers L, n c , n f .Therefore, to avoid clutter we omit the subscripts on the right hand side of eq.(10).

C. Renormalization
We regularize ultraviolet (UV) and infrared (IR) divergences of loop amplitudes in the 't Hooft-Veltman scheme of dimensional regularization, setting the space-time dimensions to D = 4 − 2ϵ.The UV divergences are removed by renormalization of the bare QCD coupling in the MS scheme.To achieve this we perform the following replacement in eq. ( 9), where S ϵ = (4π) ϵ e −ϵγ E , with γ E = −Γ ′ (1) the Euler-Mascheroni constant, and µ 0 , µ are regularization and renormalization scale respectively.The QCD β-function coefficients are The renormalized amplitudes are expanded through the renormalized coupling as in eq. ( 9).The remaining infrared divergences can be extracted through the universal factorization [48][49][50][51]: where the finite remainder R is obtained through the application of the color-space operator Z.The latter is obtained [50] from the path-ordered evolution of the anomalous dimension matrix where the sum runs over all pairs of external gluons, and the color operators T i act on the color representation of the i th parton.For adjoint indices the action is given by (T a i ) bc = −if abc .The anomalous dimensions γ cusp and γ g can be found in [52, Appendix A] 2 .
After UV and IR renormalization of amplitudes through eqs.(11) and (13) we recover expansions of eqs.( 5), ( 9) and (10) for the finite remainders R, and therefore obtain partial finite remainders which are the elementary building blocks that we focus on in this work.It is worth noting that the finite remainders contain the complete physical information about the underlying scattering process.In particular, any observable can be calculated through finite remainders (see e.g.[53]), which allows one to cancel much of the disruptions caused by the use of dimensional regularization.Finite remainders in a different IR renormalization scheme (with a different operator Z) can be obtained by an additional finite renormalization after eq. ( 13).We elaborate on this in the forthcoming publication [10].

D. Generating set of finite remainders
To calculate arbitrary observable quantities we must know all partial remainders from eq. ( 10) in all permutations in eq. ( 5), with 2 5 helicity assignments for each of them.Fortunately, the combinatorial complexity can be sidestepped by mapping each partial helicity remainders onto a small generating set.
Next, we work out the generating set for each partial remainder constructively by starting with the set of all helicities and permutations, and partitioning it into the orbits under the action of the group P ⊗ C ⊗ Σ i , where P and C are parity-and chargeconjugation respectively, and Σ i is the symmetry group of the corresponding color structure in eq. ( 5).We can then pick one representative from each orbit, and obtain the complete set of remainders by relabeling momenta and symmetries.Let us note, that for the reasons that will become clear in section III D, we do not choose the identity permutation {1, 2, 3, 4, 5} for all our representatives.Instead we prioritize to have uniform spinor weight for each of the three helicity assignments + + + + + (all-plus), + + + + − (single-minus), and + + + − − (MHV).
For the all-plus helicity configuration we have the generating remainders For the single-minus helicity configuration we have Finally, for the MHV configurations we have five generating remainders, Here we suppress the labels L, n c , n f , and, for better readability, we show the helicity labels as superscripts of the momentum labels.

E. NNLO hard function
To calculate the double-virtual contribution to NNLO QCD partonic cross sections one must square the helicity finite remainder R ⃗ h,⃗ a in eq. ( 5) and perform summation over color and helicity indices.
At leading order in α s we define the function We then define the hard function H as which can be expanded perturbatively up to O α 2 s as in eq. ( 9).Similar to eq. ( 10), we can expand H in powers of N c and N f .Through two loops we get 1) , Here again only the functions H (L),(nc,n f ) with L = n c + n f contribute in the leading-color approximation.
To perform the summations in eq. ( 21) we use the maps inverse to the ones that were used in the previous section to construct the generating set of finite remainders.It is worth noting that the polynomial expansion in N c and N f given in eq. ( 22) holds only when identities between partial remainders are correctly taken into account.

III. ANALYTIC RECONSTRUCTION
We now discuss the computation of the finite remainders (13).They can be represented as a linear combination of transcendental integral functions h i and rational coefficient functions r i , For the integral functions h i , we employ the set of non-planar pentagon functions from ref. [15].The rational coefficient functions are the central result of this paper and we obtain them via analytic reconstruction, i.e. we start with an ansatz for the rational coefficient functions r i and determine its parameters from exact numerical evaluations of the remainder over prime fields.The reconstruction cost is dominated by the time of sampling the remainders.This motivates us to search for a strategy to constrain the ansatz using physics arguments and reduce its free parameters.We will now discuss the details of this computation.

A. Numerical sampling
For numerical amplitude evaluations in a finite field we use the program Caravel [34], which implements the multi-loop numerical unitarity method [31][32][33].In this approach amplitudes are reduced to a set of master integrals by matching numerical evaluations of generalized unitarity cuts to a parametrization of the loop integrands.For the five-gluon process we use the recently obtained parametrization [35], which we extend in loopmomentum degree to match the corresponding dependence of the gluon cut diagrams.Furthermore, we extended the set of planar unitarity cuts to nonplanar diagrams which are required for subleadingcolor partial amplitudes.We generated the cut diagrams with qgraf [56] and arranged them into a hierarchy of cuts with a private code.We matched the cuts evaluated through color-ordered tree amplitudes to the amplitude definitions in section II B, by employing the unitarity based color decomposition [57,58].We extracted the ϵ-dependence of cuts that originates from the state sums in loops through the dimensional reduction method [59,60].
With these extensions Caravel now computes the integral coefficients r i of five-gluon partial amplitudes up to two loops, given a kinematic point and a choice of polarization labels for the external gluons.
Analytic expressions for the coefficient functions r i can then be reconstructed using multivariate functional reconstruction techniques [16,17] (see also recent refs.[18][19][20]25]) based on Newton and Thiele's interpolation algorithms.However, our approach is in fact closer to the ansatz-based approach of refs.[21,26].This approach constructs an ansatz for the rational integral coefficients r i , which is constrained by information from the neighborhood of their singularities.Nevertheless, here we differ even from this approach in that we use information about the residues to build linear transformations of rational functions r i to bases of functions ri with a simplified pole structure and fewer ansatz parameters.This basis change is determined numerically and simplifies the subsequent ansatz construction and parameter determination.
We will require two types of numerical evaluations for the remainder functions: 1. Random phase-space points: these are N randomly generated phase-space points which we label by the superscript n.We represent these points in terms of sets of spinor variables, They are subject to momentum conservation 2. A family of (anti-)holomorphic slices: these are Ñ < N holomorphic slices [23,35,61] associated to a subset of the random phase-space points (24), Here the reference spinor η ñ is chosen randomly.The label ñ runs over Ñ values.Similarly we will use anti-holomorphic slices which are obtained from (25) by swapping λ ñ i ↔ λñ i and η ñ → ηñ .

B. Coefficient-function basis
We now identify a basis of coefficient functions r i based on a measure for the functions' complexity.We start from the general form of the rational coefficient functions in spinor variables, The denominator factors D j are given by the symbol alphabet with integer exponents q ij [24].Computing r i amounts to determining the q ij and all parameters in the numerator polynomial N i .While it is straightforward to determine q ij (see below) it is non-trivial to determine N i because of the typically large polynomial degree.Consequently, a useful measure of complexity is the mass dimension of N i , which is linked to the number of parameters the polynomial depends on.
Simple dimensional analysis allows us to determine the mass dimension and spinor weights of the numerators N i , from those of the polynomials D j , the data q ij and the overall mass dimension and spinor weights of the helicity remainder.Consequently, we only need to determine the exponents q ij , which we compute following the univariate-slice reconstruction [24] in spinor-helicity variables [61] (see also refs.[23,35]).Here the remainder functions are reconstructed on a single holomorphic and a single anti-holomorphic slice (25), and subsequently the denominators are matched to products of the letter polynomials D j (t).This uniquely fixes the exponents q ij for each function r i .
For the five-gluon finite remainders we observe that the set of denominator factors contains the 35 elements, where the set runs over all independent permutations of the spinor strings/brackets.We also note that no coefficient in the finite remainder has a tr 5 singularity.From the q ij we can immediately deduce the mass dimension of the numerators N i .Our next goal is to determine a minimal function basis for the set of coefficient functions r i in order to reduce the number of rational functions r i that must to be constructed [60].To this end, we exploit the freedom in our choice of basis functions and select a basis containing numerators with the lowest mass dimensions.We express the linear dependent r i through the basis as where M ij is a constant rectangular matrix.At this point we do not have analytic expressions for the r i available and we require a numerical method to identify linear dependence of functions.This is achieved using a set of random evaluations (24) [60,62], which allows to associate vectors ⃗ r i of function values to each r i .A function basis is then identified by identifying a basis in the vector space of function values ⃗ r i = j∈basis ⃗ r j M ji by linear algebra.After this step, we arrive at the form of the remainder, Obtaining the analytic form of the remainder now amounts to computing the set {r i } i∈basis .

C. Basis change
So far we have selected a convenient basis of functions {r i } i∈basis .Next, we will exploit universal properties of the functions' poles, namely correlations between residues, to construct linear transformations to a simpler basis of functions, which we will denote by ri .('Simpler' again refers to the ri functions having numerators with lower mass dimension than the ones of the r i .)A basis change of this type was used already some time ago in ref. [24] and a detailed discussion of the algorithm will be presented in the forthcoming paper [36].Here we summarize the main steps.The reason that such linear combinations exist, stems from the fact that many poles (at zeros of the D i ) are spurious and cancel in the remainders.(Examples of spurious-pole denominators are spinor strings ⟨i|j + k|l], as well as higher-order spinor products [ij] and ⟨ij⟩ which do not contribute to factorization poles.)Such cancellations require that the residues of distinct r i are linearly related; only in this way they may cancel when multiplied with (degenerate) transcendental functions h i evaluated on the respective singular surfaces.
We now discuss the simplifying basis change from the functions r i to ri for a given remainder.We start by defining our requirement for a constant linear ba- Our objective is to ensure that the mass dimension of all new numerators Ñi , ri = Ñi (λ, λ) is lower than those in the original least common denominator (LCD) form (26), where dim(N i ) denotes the mass dimension of numerator N i .
The property we use is that the mass dimension is linked to the exponents qij .We are thus lead to consider the residues of the coefficient functions on the zeros of the denominator factors D j .To this end we use a family of holomorphic univariate slices (25) and obtain a univariate representation of the coefficient functions ri (t) .Near a single zero t k of one of the denominator factors we obtain a formal Laurent series with e k im being functions of the external kinematics.The e k im will be referred to as codimension-one residues.To write the above equation (34) in a more uniform way, we introduce the maximal pole degree of a given factor D k , At the same time we introduce vanishing residues, and rewrite (34) as, In this way the summation is independent of the residue function and the degree of the pole is encoded in the vanishing of residues.Equipped with this notation, we now return to the discussion of the constant transformation matrix O ij .We will now constrain the matrix, such that some of the leading residues vanish for the new basis ri .This is equivalent to lowering the respective pole degrees qij (and consequently the numerators' mass dimensions), which is what we aimed to achieve in the first place.The Laurent expansion of the functions ri depends linearly on the data of the Laurent expansion of r i , ri (t) = Consequently, we can choose the basis change O ij in such a way, that the leading residue vanishes, effectively reducing the power of the leading pole, We find that the requirement for a good basis change is that the rows of the matrix O ij are in the kernel of the leading residue functions.Independently, we also have to ensure that the basis change is invertible.Given this property, the pole degrees of the ri are reduced.Finally, we adjust the criterion for the fact that we do not have analytic expressions for the residue functions available.We follow the strategy of ref. [60,62] (see also section III B).We repeat the above construction for D holomorphic slices (25) (labeled by ñ) and obtain a representation of the residue functions e k im as vectors of their function values The linear constraint for the basis change and relies solely on numerical input.The required number D of slices is given by the dimension of the largest vector space of codimension-one residues e k jq k .Each row of O ij lives in the nullspace of ⃗ e k jq k .So far we were concerned with removing the leading singularity associated to the pole labeled by k.In reality, we impose the linear constraint (42) for multiple residues simultaneously.This includes originally subleading residues, if they transform into leading residues following the cancellation of the previous leading ones.In fact, for the most complicated function, k and m run over up to 40 residues and orders.Furthermore, we find it convenient to derive the basis change by considering one function ri at a time, which amounts to constructing O ij row by row.We start from the simplest functions, progressing towards the most complex ones.
In practice, to build each row, we perform an exhaustive breadth-first search in the space of intersections of nullspaces of the ⃗ e k im .That is, we build ri by attempting to remove an increasingly higher number of poles, until no combination exists that drops more.Figure 2 shows a simplified example of such a search tree, in the case where the function has only two denominator factors.The tree is truncated at depth 2. A tick (✓) or cross (✗) on the edge labels whether a linear combination of the r i 's exists that removes the given denominator factor.We label the numerators by their respective denominator factor powers.We look for a global minimum in the search tree.Although the D i may differ in mass dimension, we find it most convenient to minimize the sum of their exponents qij .
To ensure the rank stays maximal, after each row of O ij has been computed, the following search for ri+1 is done using ri ′ ≤i and r i ′ >i .A suitable pivoting strategy in the computation of the nullspace intersections ensures ri+1 remains a pivot, and hence the rank does not drop.Otherwise, the linear combination is discarded.In certain cases, the computation for the last rows of O ij was unnecessary, as these functions were derived through symmetries (as discussed below) or other partial remainders.The impact of the basis change is presented in table II in terms of LCD ansatz sizes.
Using basis changes that decorrelate the functions poles we have achieved much simplified coefficient functions prior to performing a prohibitively expensive reconstruction.The effect is somewhat similar to what was achieved by univariate partial fractioning by a suitable variable [22,23,63].However, we avoid introducing new spurious denominators, and in addition the order of unphysical poles is systematically reduced.Both effects are expected to lead to much improved numerical stability.Moreover, we also note that, following the change in the set of basis functions, further simplification may be achieved by a partial-fraction decomposition, either by the semi-numerical slicing techniques just mentioned, or through the purely numerical approach of refs.[21,26].

D. Analytic reconstruction and simplification
We are now in a position to reconstruct the coefficient functions.We do this by constructing an ansatz matching the pole structure of the final LCD form (31). Given previous experience, we further simplify the ansatz by removing terms with more than one ⟨i|j + k|i]-factors in the denominator [35].For example, we expect (and verify) N to be such that the following type of equality holds irrespective of what other denominator factors appear and of the degree of the spinor chains.This allows us to make the ansatz in terms of the two lowerdegree numerator polynomials N a,b instead of N .Constructing an ansatz for the numerator polynomial can be non trivial.We build ansätze for the numerator polynomials in terms of independent monomials of spinor brackets, which have both the right mass dimension and little group weights.Mathematically, this amounts to the enumeration of members of a polynomial quotient ring, subject to irreducibility by a Gröbner basis and degree bounds [21, section 2.2].The ansatz construction relies on the opensource programs Singular [64], for the Gröbner basis computation, and OR-tools [65], for the linear programming.Finally, we determine the numerator polynomials by solving linear systems for a sufficient number of random evaluations (24).
When reconstructing the coefficient functions of all remainders of the generating sets ( 17), ( 18) and ( 19), we find it most effective to reconstruct remainder functions in bunches of the same spinor weight.In fact, the permutations of momenta chosen within the above 3 sets are such that all remainders within a given set have the same spinor weight.This is useful because there exists overlap between the vector spaces of rational functions of each separate remainder.Namely, the dimension of the sum of two such vector spaces is often smaller than the sum of their dimensions.Furthermore, we find it most effective to reconstruct denominators of lower mass dimension first.We also observe that permuting the momentum assignments of the functions ri that yield the same spinor weight often yields linearly independent coefficient functions.We thus, after each newly reconstructed coefficient function, ensure that the set of functions is closed under all momentum permutations that leave the little-group weights unchanged.This ensures that information is recycled among different partials, and that the symmetries of amplitudes are exploited.Numerically, the closure of the space of functions is conveniently checked by adding all permutations and evaluating them on the random phase-space points (24).We then filter out redundant ones via Gaussian elimination.
After the reconstruction of the function coefficients is complete, we simplify the final results considering each of the three helicity sets individually.In fact, we perform a further basis change on the vector space obtained from the union of all remainders.We determine a simple basis for this space, by placing the chosen basis functions in global, or failing that local, minima in the space of possible denominator powers [35].Finally, the basis functions are then partial fractioned following the approach of refs.[21,26].
We display the set of all coefficient functions in appendices A to C.

E. Implementation
To perform the reconstruction, for the most complicated remainder, we use roughly 35,000 numerical samples.Most of these samples are points on slices needed to perform the basis change, while around 5,000 are random phase-space points.All rational functions are fitted with a single finite field.This is possible thanks to the size of the rational numbers appearing in the basis, with the largest ones not exceeding 2 digits.Initially the matrices M ij are obtained in a finite field.To lift the matrices M ij from a finite field to rational numbers we require, on the first finite field value, as many samples as the dimension of the vector space, while subsequent finite field values requires roughly a factor of 5 fewer samples for each iteration.In the end, a single evaluation with a finite field value not employed in the reconstruction is used as a check.
Finally, we summarize software packages used in the analytic reconstruction.For convenience, we take advantage of hardware acceleration on NVIDIA GPUs with the private code linac (LINear Algebra with Cuda), which we use to solve linear systems over a finite fields for the ansatz coefficients, and to handle the vector spaces of rational functions.Nevertheless, given the relatively small size of the systems, a low-level CPU implementation may also have sufficed.We use lips and pyadic [66] for the generation and manipulations of phase-space points defined in terms of spinors, and for numerical eval- Table I: For each helicity configuration, this table shows the dimension of the vector space of rational functions, and the number of functions in the generating set that spans the space upon closure under the symmetries of the helicity vector.

IV. RESULTS
We express remainders in terms of three rationalfunction bases.The function bases are obtained from a generating set of spinor-helicity functions and symmetry operations.We denote the generating functions as rh4h5 i , with h 4 and h 5 labeling the three helicity configuration: all-plus (r ++ i ), single-minus (r +− i ), and MHV (r −− i ).The respective basis functions are given in appendices A to C. The functions themselves are expressed in terms of spinor-helicity variables, and symmetry operations.The latter take the form where "−" denotes that the expression with the permuted labels should be subtracted.
Within a rational function, the employed convention is to apply the (anti-)symmetrization to all terms preceding the mapping [26, section 4.3].
To obtain the function basis from the generating functions, the set needs to be closed under the symmetries of the little group weights.These are S 5 (1,2,3,4,5), S 4 (1, 2, 3, 4) and S 3 (1, 2, 3)⊗S 2 (4,5), for all-plus, single-minus and MHV configurations, respectively (here S n (1, ..., n) denotes the group of permutations of the elements {1, ..., n}).Table I shows the dimensions of the three vector spaces, and the number of generating spinor-helicity expressions in the chosen basis.The dimensions of the vector spaces are uniquely defined, in the sense that they are the smallest spaces that are both closed under the symmetries and that span all coefficients in the Table II: For each partial amplitude, this table shows the dimension of the vector space of rational functions, and ansatz size in LCD form of the most complicated function in the basis, before and after basis change.In the cases denoted by N/A the basis change was not required.
required partials.The number of generating functions is representation dependent.
We note that the overlap between the rationalfunction spaces of the different partials is significant.This can be observed by comparing the sizes of the vector spaces in table II, to that of their union, closed under S 3 (1, 2, 3) ⊗ S 2 (4, 5), in table I, in spite of the fact that the spaces of partials are not closed under all symmetries of the phase weights.
The remaining helicity configurations are obtained by re-assigning momentum labels and/or parity conjugation.However, some of the permutation of eq. ( 5) will involve exchanges between momenta in the initial and final state.This encompasses nontrivial analytic continuation, which we perform following ref.[5].More explicitly, the action of permutation σ on a finite remainder R is given by σ Here we rely on the closure of the pentagon functions under permutations [15].The analytic continuation therefore amounts to simply obtaining a matrix M ′ ij for each required permutation.In practice, these matrices are obtained by permuting the legs in the master integrals only, and then re-mapping them to pentagon functions.

A. Ancillary files
We provide expressions for the finite remainders of all independent partial amplitudes through two loops in the decomposition of eqs.( 5), ( 9) and (10) in [67].For each helicity configuration, all-plus, singleminus and MHV, we organize the results in terms of two global bases, valid for all partials and crossings, in the files 1. basis transcendental , 2. basis rational .
The constant matrices M ij of rational numbers are partial remainder, and permutation specific.They are organized in subfolders labeling the partial remainders from subsection II D, with the notation where ⃗ h, L, n c and n f refer to helicities, number of loops, number of N c powers and number of N f powers, as defined earlier in the paper.For completeness of the ancillaries, we also provide information about permutations and color structures of the partials in the file amp info.
The matrices themselves are stored in the files rational matrix {permutation} , where it is understood that this permutation has to be combined with that of each partial as defined in subsection II D and in amp info.The order of the permutations matters, with those defined in II D taking precedence.

B. Validation
Our computation incorporates multiple internal consistency checks at various stages.In constructing the finite remainders, we ensure the cancellation of poles in the dimensional regulator ϵ at every kinematic point.Further validation of the finite field reconstruction occurs at an independent kinematic point, not utilized for the reconstruction, and with a distinct value of the prime.
We found agreement with the numerical evaluations for the helicity-and color-summed squared remainders in the leading-color approximation [5].Furthermore, we conduct additional checks through independent computations in full color.We compared the one-loop hard functions against numerical evaluations by BlackHat [68] and found agreement.At two loops we performed verification against the full-color all-plus calculation [37], and the validation of the evaluations provided in appendix D with an independent computation [40,69].

C. Numerical evaluation
We implement our analytic results for the partial helicity remainders, as well as for the NNLO hard function defined in eq. ( 21) in the C++ library FivePointAmplitudes [70], which employs PentagonFunctions++ [15,71,72] for numerical evaluation of the transcendental integral functions.This allows us to ensure the stability of numerical evaluations via the rescue system developed in ref. [5].We achieve excellent numerical performance, with a single numerical evaluation of the two-loop hard function taking a few seconds in double precision on a personal desktop computer.The average evaluation time per phase-space point on cluster nodes over the phase-space of ref. [5] is around 8s with the rescue system enabled.For comparison, the average evaluation time for the leading-color contribution only is around 1.6s with the same setup.

V. CONCLUSIONS
In this work we computed the two-loop five-gluon helicity amplitudes in QCD.The amplitudes are uniformly represented in terms of a basis of rational coefficient functions, transcendental functions and constant matrices of rational numbers that link the two.We achieved an unprecedented level of simplicity in the basis of coefficient functions, by identifying a distinguished basis of functions based on their singularity structure.However, it has to be noted that significant complexity still remains in the transcendental functions and in the mentioned matrices of rational numbers.This observation raises an intriguing possibility of a more physically motivated basis of transcendental functions, rooted in the cancellation of spurious singularities, that is, in the amplitudes' locality.Further exploration into this avenue remains an intriguing direction for future research.
In terms of phenomenological applications, the significance of our result becomes particularly evident when considering not just three-jet production at NNLO, but also two-jet at N 3 LO in hadron collisions.In fact, the latter will require integrating the provided expressions in unresolved phase-space configurations.We expect our compact basis of rational coefficient functions to benefit both precision and stability of such computations.
Finally, the simplicity achieved in the present computation is not confined to the specific channel under consideration, nor are the techniques employed specific to five-point massless kinematics.Analogous results for the quark channels of threejet production will follow shortly in a separate publication, and preliminary studies show similar benefits in tackling processes involving more challenging kinematics, such as five-point one-mass.In summary, we believe that the form of the rational coefficient functions, which directly benefits phenomenological studies, deserves further studies concerning the mathematical structure of scattering amplitudes, and may open new paths to precision predictions for multi-scale scattering processes at particle colliders.

ACKNOWLEDGMENTS
First and foremost, we gratefully acknowledge Ben Page for stimulating conversations and for collaboration on the basis-change algorithm [36] that we applied in this paper.We thank Thomas Gehrmann for interesting discussions.We thank Johannes Schlenk for discussions and comments on the manuscript.Finally, we thank Bakul Agarwal, Federico Buccioni, Federica Devoto, Giulio Gambuti, Andreas von Manteuffel, Lorenzo Tancredi for the numerical comparison of the reference values which are displayed in appendix D. We gratefully acknowledge the computing resources provided by the Paul Scherrer Insitut (PSI) and the University of Zurich (UZH).V.S. has received funding from the European Re-search Council (ERC) under the European Union's Horizon 2020 research and innovation programme grant agreement 101019620 (ERC Advanced Grant TOPUP).G.D.L.'s work is supported in part by the U.K. Royal Society through Grant URF\R1\20109.

Figure 2 :
Figure 2: Example of a simple search tree for intersecting null spaces with a global (r i ) and a local minimum ( Ñ2,1 /(D 2 1 D 2 )).