Two-loop leading colour QCD corrections to $q \bar{q} \to \gamma \gamma g$ and $q g \to \gamma \gamma q$

We present the leading colour and light fermionic planar two-loop corrections for the production of two photons and a jet in the quark-antiquark and quark-gluon channels. In particular, we compute the interference of the two-loop amplitudes with the corresponding tree level ones, summed over colours and polarisations. Our calculation uses the latest advancements in the algorithms for integration-by-parts reduction and multivariate partial fraction decomposition to produce compact and easy-to-use results. We have implemented our results in an efficient C++ numerical code. We also provide their analytic expressions in Mathematica format.


Introduction
The production of a pair of isolated highly-energetic photons in proton-proton collisions, pp → γγ + X, represents an important class of processes for the physics programme at the Large Hadron Collider (LHC) and at future hadron colliders.Among many relevant aspects, pairs of prompt photons (diphotons) constitute an irreducible background to various Standard Model (SM) and Beyond Standard Model (BSM) processes, most prominently Higgs boson production in its H → γγ decay channel.Indeed, with the increase of collision energy, the diphoton invariant mass distribution can provide a powerful tool to search for heavy resonances decaying to pairs of photons [1,2], while its transverse momentum distribution offers a unique probe to investigate their properties.It is therefore crucial to provide an accurate theoretical description of the production of a pair of photons recoiling against hard Quantum Chromodynamics (QCD) radiation, across a vast spectrum of energies.
At the technical level, the theoretical description of the production of a pair of photons with large transverse momentum is non trivial for several reasons.An obvious one is that it requires computation of 2 → n scattering amplitudes with n ≥ 3 for various partonic channels relevant at a specific perturbative order.In particular, at leading order (LO) the process receives contribution from three partonic sub-channels, q q → gγγ, qg → qγγ, and qg → qγγ, where the second and third channels can be obtained as crossings of the first.The loop-induced gluon-fusion process gg → gγγ, instead, contributes formally only at next-to-next-to-leading-order (NNLO).
The mathematical complexity of loop corrections to the scattering amplitudes above is the main reason why pp → γγ + jet is currently known only up to next-to-leading-order (NLO) QCD [3,4].Calculation of cross section and differential distributions for pp → γγ+jet through NNLO QCD requires the computation of two-loop amplitudes for q q → gγγ and its crossings, together with an efficient subtraction scheme to organise and cancel the infrared (IR) divergences between real and virtual contributions.In recent years, a lot of progress has been achieved on both fronts.On the IR subtraction side, several schemes have been developed which are able to handle NNLO QCD corrections to the production of a colour singlet plus one QCD jet [5][6][7][8][9][10][11].On the amplitude side, equally impressive results have been achieved.The first two-loop 2 → 3 scattering amplitudes for the production of three jets and three photons in leading-colour QCD have been computed [12][13][14][15][16][17][18][19][20][21][22][23].Together with the developments on the IR subtraction side, this has made it possible to provide the first NNLO studies for production of three photons at the LHC [24,25].
In this paper, we consider the calculation of two-loop QCD corrections to the production of a diphoton pair and a jet for the partonic channels q q → gγγ and qg → qγγ, based on a Feynman diagrammatic approach.We demonstrate how the use of new algorithms for the reduction of loop integrals and multivariate partial fraction decomposition allows us to compute these two-loop corrections in a rather straightforward way and to produce very compact and efficient numerical implementations for them.While it is common practice to compute helicity amplitudes for multi-loop processes, the case of diphoton production plus jet does not require us to keep track of the polarisations of the external photons.We therefore expect it to be sufficient for near-term phenomenological applications to consider the interference of the two-loop amplitudes with the corresponding tree-level amplitude, summed over colours and polarisations.We will show that, also in this case, very compact results can be obtained, similar to what can be achieved for comparable helicity amplitudes.Together with the recently computed three-loop QCD corrections to diphoton production [57], these amplitudes also provide an essential ingredient towards the calculation of pp → γγ in N 3 LO QCD.
The rest of the paper is organised as follows.In Section 2 we describe the processes and their kinematics.In Section 3 we illustrate the general structure of the scattering amplitudes, and in particular of the interference terms contributing to the squared amplitude.Technical aspects of the diagrammatic calculation, the integral reductions and the multivariate partial decompositions are presented in Section 4. In Section 5 we describe our renormalisation and infrared subtraction procedures which define the final form of the results.In Section 6 we discuss how we optimise our analytic results by using a minimal set of rational functions, and in Section 7 we present their implementation in a C++ numerical code.We finally draw our conclusions in Section 8.

Kinematics and Notation
We consider the production of a pair of photons in association with a gluon in quarkantiquark annihilation up to two-loop order in massless QCD, i.e. corrections up to O(α 2 s ) relative to the treelevel.We focus here on the scattering process in the physical region, i.e. a 2 → 3 kinematic configuration.For simplicity, in what follows, we do not discuss in detail the other partonic subchannel, qg → qγγ, which can be obtained by crossing symmetry according to q(p 1 ) + g(p 2 ) → q(p 3 ) + γ(p 4 ) + γ(p 5 ) . (2.2) Analytic results for both channels, with the appropriate identification of the external momenta as in Eqs.(2.1) and (2.2), are provided in the ancillary files of the arXiv submission of this paper.Kinematics are fixed by imposing that all external particles fulfil the on-shell condition p 2 i = 0, such that one is left with five independent kinematic invariants, which we choose to be the five adjacent ones, It is also useful to introduce the parity-odd invariant where µνρσ is the totally anti-symmetric Levi-Civita symbol.The invariant 5 is related to the determinant of the Gram matrix G ij through where we find it convenient to factor out the SU(3) colour generator T a ij , with i and j being the colour indices of the quark and anti-quark and a the colour index of the gluon in the adjoint representation.* (p 3 ) is the polarisation vector of the outgoing gluon and * (p 4 ), * (p 5 ) are the ones of the photons.Clearly, for the qg channel in (2.2) the same applies upon suitably renaming the external momenta p 2 ↔ p 3 and ommitting the complex conjugate for the polarisation vector of the incoming gluon.The amplitude A stripped of the colour generator T a is then perturbatively expanded in the strong coupling constant α s as with α the fine-structure constant and Q q the electric charge of the quarks in units of electron charge.The decomposition of (3.2) fixes the normalisation of the expansion terms A i .The amplitude squared and summed over colours and polarisations can be expressed as col,pol where we introduced the interference terms and a global factor C, that accounts for colour sums, which is defined as (3.5)Note that all colour and helicity degrees of freedom are summed over in (3.3) rather than performing an average for the initial states.The perturbative corrections to the expansion coefficients for the amplitudes A i , and to their interferences W ij , can be decomposed in terms of their colour and electroweak factors which can be expressed as polynomials in where C A and C F are the quadratic Casimirs for the colour group SU(N ), (N = 3 for QCD), and n f is the number of (massless) quarks running in closed fermionic loops.The factors n γ f and n γγ f correspond to the number of quarks, weighted with their electric charge, running in closed fermion loops, with one and two external photons attached to, respectively.The coefficients of various powers of the factors in (3.6) constitute gauge-invariant subsets of the final result, thus it is natural to decompose the amplitude in terms of them.
Within the prescription of (3.4)-(3.5), the LO result W 00 is a rational function of invariants only, while the tree with one-loop interference admits the decomposition The last contribution is introduced by the renormalisation procedure, see Section 5 for details.The tree-two loop interference has a richer structure, which we decide to organise as with the individual colour factors given by where c 9 arises from a contraction of type d abc d abc and c 10 follows again from ultraviolet (UV) renormalisation.A selection of representative diagrams contributing to each colour structure is shown in Fig. 1, where diagram (i) contributes to the factor c i and the colour factor c 10 follows from diagrams of type (9) as well.Eventually, one can then expand the polynomial in N for c 1,2,3 and write where W (1) 02 is the leading colour (LC), W 02 is the next-to-leading colour (NLC) and W (3) 02 is the next-to-next-to-leading colour (NNLC) contribution.The expression in terms of Casimir operators can always be recovered through the set of identities In this paper we will focus on the calculation of the LC contribution W (1) 02 and the three fermionic contributions W (i) 02 , i = 4, 5, 6, which can be obtained by considering only the planar two-loop diagrams.

Diagrammatic calculation and integral reductions
We perform our calculation in conventional dimensional regularisation (CDR) [59,60], working in d = 4 − 2 dimensions.UV and IR singularities will then manifest as poles in the regulator .
In order to compute the coefficients W 00 , W 01 and W 02 , we start by producing all relevant Feynman diagrams with QGRAF [61].We find 6 diagrams at tree level, 80 diagrams at one loop and 1716 diagrams at two loops.We use FORM [62] to manipulate the diagrams, interfere them with their tree-level counterparts and perform the relevant Lorentz and Dirac algebra.Starting at two loops, it is particularly convenient to group the diagrams depending on the graph they can be mapped to, before performing all symbolic manipulations.This allows us to avoid repeating expensive operations multiple times.After this step is complete, we can express our two-loop interfered amplitude in terms of scalar Feynman integrals drawn from two different integral families.We define the integrals as where γ E ∼ 0.5772 is the Euler constant and, following the notation of [44,45], we define the two families in Table 1.
As a first step to simplify our interference terms, we use Reduze 2 [63,64] to search for symmetry relations among the different scalar integrals in the two topologies in Table 1 and their crossings.This allows us to substantially reduce the number of different integrals that we need to reduce to master integrals.First, we collect the integrals that are required to express the coefficients of any of the colour factors (3.9).We find that after applying the symmetrisations above, and modulo crossings of the external legs, the interference of the two-loop amplitude with the tree-level amplitude requires the reduction of 1811 integrals in family A and 508 integrals in family B, which include all the integrals required to Prop.den.Family A Family B Table 1.Definition of the two integral families used in the calculation.We note that also both the non-planar hexagon-box and the double-pentagon topologies can be described by family B and crossings thereof.
compute both the hexagon-box and the double-pentagon topologies.In both families we need at most rank-5 integrals (up to 5 irreducible scalar products in the numerator).From this point on, we focus only on the planar colour factors and therefore only on integrals belonging to family A. We performed their reduction using both Kira [35,65,66] and an in-house implementation, Finred, which employs finite field sampling, syzygy techniques, and denominator guessing, see [37,67] for more details.It is worth stressing that the planar reductions up to rank 5 did not constitute an issue for either program here, and took e.g.40 hours on 36 cores with Kira and a similar runtime with Finred.The reduction lists for family A produced in this way are not extremely complicated, with a size of around 390 MB in total.Note that we performed the reduction directly in terms of the pre-canonical set of master integrals defined in [44].We stress here that these lists do not include the crossings necessary to reduce all diagrams.
To arrive at a complete set of reduction identities and to render their inclusion as simple as possible, we proceed as follows.The integration-by-parts solvers deliver each integral coefficient as a rational function in a common-denominator representation.We find it useful to convert the rational functions to a partial fraction decomposed form.Due to the choice of master integrals, we encounter only irreducible denominator factors which depend on either d or the kinematic variables.For the d dependence of the rational functions, a simple univariate partial fraction decomposition is sufficient.In contrast, the decomposition involving the kinematic denominators is computationally non-trivial.We For them, we employ the MultivariateApart package [67] to perform a multivariate partial fraction decomposition; see also [68][69][70] for related decomposition techniques.The decomposition of [67] is based on replacing all irreducible denominator factors d k ({s ij }) in a given expression according to and reducing the resulting polynomial with respect to the ideal generated by using a suitable monomial ordering.Due to the specific structure of our denominator list, the monomial block ordering of [67] coincides with a lexicographic ordering of the q i before any s ij are considered, which produces a Leȋnartas decomposition [71,72].Note that we included only the local denominators of each coefficient in the initial partial fraction decomposition.This run took a couple of days and reduced the size of the reduction list by almost a factor 5, down to around 80 MB.In particular, the reduction in complexity is particularly pronounced for the most complicated rank-5 identities, for which up to a factor of 40 reduction in size is seen.While it has been known for a long time that integration-by-parts identities can become substantially simpler even using naive variants of multivariate partial fraction decompositions, we would like to point out a systematic study of the impact of these new algorithms on the size of the reduction identities which has recently appeared in [70].Starting from these substantially simpler identities, we perform the relevant crossings and then a second partial fraction decomposition.The second partial fraction decomposition employs a global Gröbner basis for the ideal (4.4), taking into account the denominator factors of all coefficients at once.We emphasize that this method allows to decompose terms of a sum locally, but still ensures a globally unique representation of the rational functions in the results.With this multi-step procedure, we obtained all identities and their crossings in a form suitable for insertion in the diagrammatic calculation, while keeping the dimension of the expressions under control at each step.After insertion of the reduction identities into the interference terms, a final quick last partial fraction decomposition handles the remaining few additional denominator factors.In practice, we find the partial fraction decomposition of the original, uncrossed identities to be the most time-consuming step, which nevertheless could be handled in a couple of days.The remaining steps to arrive at the reduced and fully partial fractioned interference terms took a few hours.
In order to optimize our results for numerical stability, we found it useful to tune our monomial ordering for the partial fraction decomposition in such a way that spurious poles and unnecessarily high powers of the denominators are avoided.This choice of ordering was performed for each colour factor and loop order separately.In this way, we obtained a rather compact expression for the two-loop interference terms in terms of canonical master integrals, with full dependence on the space-time dimensions d.As we will see in Section 6, by substituting the explicit results for the master integrals, expanding about d = 4, subtracting the poles, and properly simplifying the remaining rational functions, we are able to obtain extremely compact, fully analytic results for the finite remainders of the interference terms.

Renormalisation and infrared factorisation
We work in fully massless QCD and keep the number of light fermion flavours n f generic; we do not consider any loop corrections due to massive quarks.We perform UV renormalisation in the standard MS scheme, which allows one to express renormalised amplitudes in terms of unrenormalised ones by simply replacing the bare coupling constant α b s with the running coupling α s (µ 2 ), evaluated at the scale µ 2 , where S = (4π) e − γ E , and β 0 and β 1 denote the first two perturbative orders of the QCD beta function, with T r = 1/2.The renormalised interference terms in (3.3) are explicitly obtained as where with a superscript b we denote bare quantities.After UV renormalisation the results in (5.3) contain only IR poles.We subtract them according to Catani's factorisation formula [73], which makes it possible to reorganise the interference terms as where the action of the operators I 1 and I 2 produce the complete IR structure of the renormalised amplitudes.The finite remainders W fin 01 , W fin 02 will be the main result of our paper.Before discussing these in detail, let us show the explicit formulae for the Catani operators for the process under consideration. 1 By direct calculation it is easy to see that for the q q channel where γ g = β 0 , and where K is universal, and H 2 is a renormalisation and process-dependent factor, in our case given by with H 2,q and H 2,g which read explicitly [74] H (5.10) In order to obtain the corresponding expressions for the qg channel, it is enough to swap the indices 2 ↔ 3 in (5.5) and the rest follows.
Since we are working in CDR, we carry out the computation of the LO term W 00 retaining exact d dependence and we expand W 01 up to O( 2 ).This is necessary for the correct subtraction of the IR poles according to (5.4).
After UV renormalisation and IR subtraction the one-loop finite remainder takes the form: W fin 01 = C A W (1),fin 01 where, as anticipated in Section 3, the last term on the right-hand side originates from the renormalisation and subtraction procedures and it is proportional to the LO result.As for the two-loop finite remainder, we find it convenient to cast it in the form where the colour factors c 4,...,10 are the same as defined in (3.9) and we introduced (5.13) The factors c1,2,3 clearly organise the colour tower, while c 10 is introduced by UV renormalisation and IR subtraction.As we have already argued at the end of Sec. 3, the coefficients W (1),fin 02 , W ,fin 02 , W (5),fin 02 and W (6),fin 02 entail only planar two-loop corrections.These four finite contributions constitute the main object of this paper.
A first test that we performed on our results at the analytical level, is that we reproduce all the poles of the renormalised interference terms, at one and two loops, as predicted by Catani's factorisation formula in (5.4).

Exploiting linear dependencies of rational functions
The LO contribution defined in (3.3) where we kept full -dependence and in the last row we denote cyclic permutations of the indices labelling the kinematic invariants.An equivalent expression holds for process (2.2) upon replacing 2 ↔ 3 in the indices labeling the invariants.
Starting already from one-loop, the explicit expression of the interference terms become too lengthy to be shown here.Therefore, we limit ourselves to describing the functional form of our results and the approach we adopted to obtain particularly compact expressions.
Both at one and two loops, the finite remainders will be a linear combination of transcendental functions multiplied by rational functions.The physical one-loop finite remainders, i.e. the coefficients of O( 0 ), contain only simple logarithms and dilogarithms and they are free of the parity odd-invariant 5 .At the two-loop order, as well as for the higher powers of the one loop result, the class of transcendental functions needs to be extended beyond classical polylogarithms and 5 enters explicitly.We employ the representation in terms of so-called Pentagon Functions presented in [47] 2 .Thus we write a generic finite remainder as where m denotes some given colour factor, n = 1, 2, r k is a rational function of the kinematic invariants and f k indicates a pentagon function.Equation (6.2) holds formally at LO, n = 0, as well, upon putting f k = 1.We remind the reader that at this step we represent each rational function r k in its partial fractions decomposed form.We further stress that the algorithm presented in [67] avoids the presence of spurious denominators, therefore our final expressions contain all and only the physical denominators of the amplitudes.Despite that the partial fraction decomposition allows us to obtain rather compact results already at this level, we find that the representation in (6.2) is not the most compact representation yet.This is due to the fact that the rational functions r k are not linearly independent and, in fact, the set of independent rational functions is substantially smaller.The minimal set of rational functions can be obtained as follows.First of all, in order to guarantee a globally unique representation of all rational functions, we employ a partial fraction decomposition with respect to a common Gröbner basis.As a result, each rational function r k is represented as a polynomial of maximum degree p in the 5 kinematic variables s ij (2.3) and the 25 inverse denominators q i (4.3), and the coefficients a k,m 1 ...m 30 are rational numbers.We emphasize that only certain M k,m 1 ...m 30 occur as irreducible monomials in our partial fractioned results; in practice it is convenient to enumerate only those.We are interested in determining linear relations between different rational functions, and employing them to re-express all r k in terms of a linearly independent subset of them.By inserting the partial fractioned form (6.3) and observing that the monomials M m 1 ...m 30 are independent, we obtain a system of linear relations for the b k , for each of the irreducible monomials M m 1 ...m 30 .In our case, there are many more monomials than rational functions, and the system is over-constrained.Nevertheless, many of the equations turn out to be linearly dependent, allowing us to find a solution to the system in terms of a basis of independent rational functions.Note that this can be achieved by a row reduction of a matrix of rational numbers, which can easily be done e.g. with a finite field solver like Finred.Similarly to what happens for the reduction to master integrals, there is not a unique solution for the basis of rational functions and a more natural choice can lead to more compact results for the scattering amplitude.We find that very compact analytic expressions can be found by simplifying the rational functions for each loop order and colour factor separately.In practice, we order the rational functions according to the number of monomials they contain, and at each step we remove the one with the largest number.This approach makes it possible to drastically reduce the size of the final results.As an example, for the LC coefficient W (1),fin 02 we start with an expression of around 64 MB which, after moving to a basis of independent rational functions, can be reduced to around 3 MB.This makes it possible not only to have a more compact result, but, most importantly, a more efficient and potentially more stable numerical implementation.Clearly, we cannot demonstrate that a more compact representation does not exist for a different choice of independent rational functions.On the contrary, we expect that more compact expressions could be obtained if one does not insist on using a set of independent kinematics invariants, see for example Eq. (6.1) and Ref. [22].Nevertheless, our current representation is more than suitable for practical use.
As last a remark, we stress that the procedure above does not produce the minimal number of rational functions.Indeed, starting from the results simplified according to the procedure above, we can attempt to further relate rational functions across different colour structures and different loop orders.In this way, more relations can be found and a minimal set of rational functions can be identified.We find this second representation particularly suitable for implementation of the interference terms in a numerical code, as described in the following section.We provide analytical expressions for various colour factors at the different loop orders with the arXiv submission of this manuscript.

Numerical implementation
We implemented our results in the C++ numerical code aajamp, which we distribute through the git repository [75]; it can be easily downloaded with git clone git@gitlab.msu.edu:vmante/aajamp.gitDetails on the installation procedure and usage of the package are given in the git repository.The main purpose of the code is to evaluate the finite remainders W 00 , all W (i),fin 01 of (5.11) and W (i),fin 02 of (5.12) for i ∈ {1, 4, 5, 6, 10}.The evaluation of all transcendental functions in our code, both at one and two-loop level, is handled by the PentagonFunctions-cpp library of [47,76,77].In order to implement our formulae in an efficient fashion we adopted the optimised code generation provided by FORM [78].
We then checked the implementation of our LO and NLO results against OpenLoops 2 [79] and find agreement.
In Table 2 we provide benchmark results for a kinematic point in the physical region defined by To assess the performance of our code, we measured the evaluation time of the NLO and NNLO results in double precision on a single Intel i7-9750H CPU @ 2.60GHz core using gcc 9.3.0 for a distribution of physical points in phase space, which we describe in more detail below.We find an average evalution time of 5.2 × 10 −2 ms and 1.2 s per phase space point for the NLO and NNLO contributions, respectively.We note that most of the evaluation time for a given point is spent on the computation of the pentagon functions.
The aajamp package allows for numerical evaluations in double and quadruple floatingpoint precision, i.e. 16 and 32 decimal digits representation.The quadruple precision evaluation is adapted to follow the strategy of the PentagonFunctions-cpp library, which  2. Benchmark results for the tree-level, one-loop and two-loop interference terms for the partonic channels q q and qg for the kinematic point in (7.1).Only the real part of W (i) 0n is presented.
utilizes the qd library [80] for higher-precision arithmetic. 3Therefore, aajamp relies on qd as well.We note that the PentagonFunctions-cpp and qd libraries offer octuple precision, but we reckon that quadruple precision is sufficient for phenomenological applications.Although the user can choose at the interface level between a purely double or quadruple arithmetic evaluation, we have implemented a primitive precision control system.This is motivated by the fact that individual terms appearing in the squared matrix elements can develop spurious singularities inside the physical phase space, leading to more or less severe numerical instabilities.These are typically associated with small Gram determinants, or in our case, with denominators in the rational functions which become small compared to the typical energy scale of the process.We observe that in our final expressions, the parity odd invariant 5 does not appear in any denominator and, correspondingly, small Gram determinants are not of prior concern.Instead, we focus on small denominators in the rational functions, which are not associated with a physical singularity, but lead to a major source of numerical instabilities in the evaluations.Let us stress that, in principle, also other types of numerical instabilities could arise, but our denominator based analysis indeed works well in practice, as we will show.Examples of these spurious singularities arise in events with all particle having a large transverse momentum and entailing collinear photon pairs or collinear gluon-photon pairs.We therefore find it natural to activate a quadruple precision evaluation if any of the 25 denominators becomes smaller than a given threshold which the user can tune.As can be seen in Fig. 2, this precision-control system significantly improves the reliability of the numerical evaluations.Here, we assess the level of numerical stability via a rescaling test, i.e. we consider the quantity D defined as .The labels "dp" and "rescue" are described in the text.Events were generated according to (7.3) where we perform two evaluations of the colour factor W (i),fin 02 with input kinematics (s ij , µ 2 ) and with the same kinematics rescaled by a factor ξ. The extra factor ξ n accounts for the mass dimensionality of W (i),fin lm , and in our case we just have n = 1.The definition of D in (7.2) intuitively provides an indicator for the number of digits one can trust for a given computation, thus we identify it as our level of instability.In Fig. 2 we show a cumulative histogram of phase-space points, which provide an evaluation of W (i),fin 02 with an instability D larger than D min .We plot two selected colour factors, both evaluated in pure double precision, labelled as "dp", and with the precision control system activated, labelled as "rescue".We generated 10 5 uniformly distributed events with Rambo [81] subject to the constraints E com = 1 TeV, p T,g > 30 GeV, p T,γ 1 > 30 GeV, p T,γ 2 > 30 GeV, (7.3) where E com is the energy in rest frame of the colliding partons and p T,i is the transverse momentum of particle i.One can see that for this ensemble of events, the precision control system is able to capture and cure the worst instabilities, effectively guaranteeing a very good level of accuracy.Once again, we stress that a more sophisticated stability system could be devised in principle, especially in regions of soft or collinear emissions.

Conclusions
In this paper, we presented the calculation of the leading colour and light fermionic twoloop corrections for the production of two photons and a jet in quark-antiquark annihilation and quark-gluon scattering.This calculation has been made possible by the combination of state-of-the-art techniques for the reduction of multiloop Feynman integrals, new ideas about their representations in terms of multivariate partial fractions, and recent results for the relevant master integrals.In particular, we have shown how the spin summed interference between the two-loop and the tree amplitude can be computed from its Feynman diagrammatic representation, resulting in very compact analytic expressions.We have checked that our two-loop corrections display the correct pole structure, as first predicted by Catani for QCD amplitudes at this perturbative order.In order to demonstrate the flexibility and usability of our result, we have also implemented the tree-level, one-loop and two-loop finite remainders in a C++ library, providing all expansions through to transcendental weight four.Our library links against the library and allows the user to evaluate the loop corrections in double and quadruple precision.We have also introduced a simple precision control system that allows the code to identify phase space points prone to loss of precision, such that quadruple precision evaluations are restricted to a minimum.We envisage that the algorithms developed for this calculation can be extended to solve future cutting-edge problems in the computation of multiloop multileg scattering amplitudes relevant for collider physics phenomenology.

Figure 1 .
Figure 1.Examples of two-loop diagrams which contribute to the various colour factors in (3.9).

Figure 2 .
Figure 2. Probability finding an event with an instability level D defined in (7.2) larger than D min for the two colour factors W (1),fin 02 and W (4),fin 02