Leading-color two-loop QCD corrections for three-photon production at hadron colliders

We compute the two-loop helicity amplitudes for the production of three photons at hadron colliders in QCD at leading-color. Using the two-loop numerical unitarity method coupled with analytic reconstruction techniques, we obtain the decomposition of the two-loop amplitudes in terms of master integrals in analytic form. These expressions are valid to all orders in the dimensional regulator. We use them to compute the two-loop finite remainders, which are given in a form that can be efficiently evaluated across the whole physical phase space. We further package these results in a public code which assembles the helicity-summed squared two-loop remainders, whose numerical stability across phase-space is demonstrated. This is the first time that a five-point two-loop process is publicly available for immediate phenomenological applications.


Introduction
Precise theoretical predictions for scattering experiments at particle colliders crucially rely on the availability of higher-order scattering amplitudes.Over the next few years, the large increase in the amount of data collected by the experiments at the Large Hadron Collider (LHC) at CERN will translate into measurements made at unprecedented accuracy.To make the most of the physics program of the LHC it is fundamental that theoretical predictions reach comparable levels of precision.A large number of 2 → 2 scattering processes are now available including next-to-next-to-leading (NNLO) order QCD corrections (see e.g.[1] for a recent review).However, reaching the same level of precision for higher multiplicity processes still remains a formidable challenge, with the best available predictions only including next-to-leading order (NLO) QCD and electro-weak corrections.
In this work, we focus on the double-virtual NNLO QCD corrections to the production of three photons at hadron colliders.This process is important for studying various beyond-the-Standard-Model (BSM) phenomena.In particular, it can be used to constrain anomalous quartic gauge [2] and Higgs couplings [3][4][5].Furthermore, it is a major contribution to the irreducible background in searches for associated production of BSM particles and a photon [6][7][8].Similar to the related two-photon production process [9][10][11][12][13][14], threephoton production exhibits very slow perturbative convergence [15,16] and non-reliable estimates of uncertainty from missing higher orders.In particular, a significant tension between LHC Run 1 data and NLO QCD predictions has been reported [17].The computation of NNLO QCD corrections is thus crucial for this process.While several NNLO infrared subtraction approaches are in principle capable of tackling 2 → 3 processes, the computation of the two-loop five-point amplitudes required for the double-virtual contributions to three-photon production remains a major obstacle.
In the last few years, there has been remarkable progress in the calculation of twoloop five-point amplitudes.A basis of master integrals relevant for the scattering of five massless particles has been known for a number of years in the planar limit [18,19], and more recently including non-planar topologies [20,21].These master integrals evaluate to a class of multi-valued special functions with logarithmic branch cuts.Given the large number of scales in five-point kinematics, it is essential to find a representation of the master integrals in terms of special functions that manifests the physical properties of amplitudes.This was achieved by defining a set of so-called pentagon functions, first for the planar integral topologies [22], and more recently for a complete set of five-point massless master integrals [23].The latter work also provides an efficient code for their numerical evaluation across all of phase space.At the same time, there has been substantial progress in the reduction of five-point two-loop scattering amplitudes to master integrals [24][25][26][27][28][29][30][31][32][33][34][35][36], often building on the use of finite-field arithmetic [26,37], and functional reconstruction techniques [26,31,[38][39][40].As a result, all five-parton two-loop planar QCD amplitudes were computed numerically [41][42][43][44][45] and in analytic form [46][47][48][49].Despite the remarkable progress of the last few years, five-point amplitudes are not yet broadly available in a form suitable for direct phenomenological applications.Indeed, the only phenomenological study involving five-point two-loop amplitudes is the work of [15], where the authors evaluated planar two-loop q q → γγγ amplitudes on a sufficient set of phase-space points to construct interpolating functions.
In this paper, we calculate a complete set of independent planar two-loop helicity amplitudes required for the double-virtual NNLO QCD corrections to three-photon production at hadron colliders.We present them for the first time in a form which is suitable for direct and flexible phenomenological applications.Indeed, our results have already been employed in ref. [16] to implement the first on-the-fly computation of NNLO QCD corrections to tri-photon production.We obtain the two-loop helicity amplitudes by following a similar approach to that used in refs.[47,49], based on the two-loop numerical unitarity approach [43,50] as implemented in the recently released C++ library Caravel [36].Unlike in refs.[47,49], where only the finite remainders were obtained, in this paper we obtain analytic expressions for the decomposition of the amplitudes into master integrals, marking the first time that such results are available for five-point two-loop amplitudes.This decomposition is valid to all orders in the dimensional regulator.They are valuable for studying the analytic complexity of the master-integral decomposition and for future computations of higher-order corrections.We also obtain analytic expressions for the twoloop finite remainders, decomposed into the pentagon functions of ref. [23].We confirm that the latter are more suitable for efficient and stable numerical evaluations.All the analytic results we obtain are made available in a set of ancillary files, and we provide a public C++ library [51] for the efficient numerical evaluation of the finite remainders.
The paper is organized as follows.In section 2 we establish our conventions and define the objects we will be computing.In section 3, we present our approach to the calculation of the two-loop amplitudes and their finite remainders.We discuss some properties of the analytic results we obtained, give reference evaluation values and discuss the checks that were performed.In section 4, we present a public implementation of the numerical evaluation of the analytic results for finite remainders and demonstrate its efficiency and numerical stability.

Helicity Amplitudes
We consider a complete set of helicity amplitudes required for the computation of the double-virtual next-to-next-to leading order (NNLO) QCD corrections to the production of three photons at hadron colliders.More precisely, we consider the parton-level scattering process This is the only (sub-)process required for these corrections, as the loop-induced process gg → γγγ vanishes to all orders in the coupling constants due to the charge conjugation symmetry of QCD⊗QED.We denote the momenta and the helicity states of the particles as p i and h i respectively.All particles are massless, and the kinematics of the process are thus fully specified by the five independent Mandelstam invariants and the parity-odd contraction of four momenta where ε( The latter condition is equivalent to the negativity of the five-point Gram determinant and it is trivially satisfied by real-valued momenta in eq.(2.3).The amplitudes for this process, denoted as M(1 h 1 q , 2 h 2 q , 3 h 3 γ , 4 h 4 γ , 5 h 5 γ ), can be decomposed into a color factor, a helicity-dependent spinor weight, and a Lorentz invariant kinematic factor.That is, we can write where i 1 and i 2 are color indices of the quark and the antiquark, e q is the electric charge of the external quark in the process, and Φ denotes the spinor-weight factor.In the following, we will call A(1 h 1 q , 2 h 2 q , 3 h 3 γ , 4 h 4 γ , 5 h 5 γ ) the helicity amplitudes for the process in eq.(2.1).For simplicity, we will often suppress the arguments of A. We employ the 't Hooft-Veltman scheme of dimensional regularization with D = 4−2 space-time dimensions to regularize infrared and ultraviolet divergences.We define the dimensionally regularized helicity amplitudes with external quarks as in [45].
The bare helicity amplitudes have a perturbative expansion in powers of the bare strong coupling α 0 s = (g 0 s ) 2 /(4π), which we write as The renormalized coupling α s is related to the bare α 0 s through where γ E is the Euler-Mascheroni constant, and µ 0 and µ are the dimensional regularization and renormalization scales, which from now on we assume to be equal.β 0 is the first coefficient of the QCD β-function, where C A = N c is the quadratic Casimir of the adjoint representation of the SU (N c ) group, and T F = 1/2 is the normalization of fundamental representation generators.Below we will also need the quadratic Casimir of the fundamental representation, C F = Nc 2 −1 2Nc .We define the perturbative expansion of the renormalized amplitudes as R + O(α 3 s ) . (2.9) The coefficients R are then related to their bare counterparts as 1) . (2.10) There are 16 different helicity configurations to consider.However, it can be easily shown that only 2 of them are independent.In this paper we choose the independent configurations to be where we indexed each independent amplitude by the photon helicities.In (2.12) Each helicity amplitude can be further decomposed into individually gauge-invariant contributions as where N f is the number of light quarks, and Q f is the ratio of the electric charges of the quark with flavor f and the quark in the initial state.Here we do not consider the contributions from heavy quark loops.In fig. 1 we depict representative diagrams for each of these contributions. 1n this paper, we compute the two-loop amplitudes A (2) in the leading-color approximation, where the number of colors N c is large with the ratio N f /N c kept constant.In this approximation only planar topologies contribute.We also do not consider the gaugeinvariant term Ã(2,N f ) which includes non-planar contributions.The validity of this approximation for phenomenology is discussed in ref. [16].We write the two-loop amplitudes A (2)  as In summary, the main goal of this paper consists in the calculation of four functions: −++ .We note that since A (0) +++ = 0, only the latter two are required to compute NNLO QCD corrections to the process in eq.(2.1).
Representative diagrams for each of the contributions in eq.(2.13).Photons (γ) are denoted by wavy lines, gluon by curly lines, and quarks by straight lines.

Finite Remainders
At NNLO, the two-loop amplitudes contribute to the cross section only through finite remainders (see e.g.[52]), which can be expanded similarly to eq. (2.9) as The coefficients R (i) are obtained from the expansion of the renormalized amplitudes A R by subtracting the remaining infrared singularities [53].More precisely, we define (2.16) where the functions I (1) and I (2) are defined as (see e.g.[9,10]) (2.17) In I (2) ( ), we have introduced the functions (2.18) The finite remainders can be decomposed similarly to eq. (2.14), and we write Analogously to eq. (2.11), whenever convenient we will also write R (i,j) +++ and R (i,j) −++ to denote the contributions to each of the two independent helicity states.

Calculation of Helicity Amplitudes
To compute the functions A −++ defined in eq.(2.14), we use the framework of two-loop numerical unitarity [25,43,50,54] coupled with functional reconstruction techniques.The same approach was already used previously to compute the five-parton two-loop amplitudes [47,49].We build on the implementation of this framework in Caravel [36], which we modify to handle amplitudes with external photons.In this section, we summarize the main steps of our calculation.
Before delving into this, however, we note that our approach requires knowing the one-loop amplitudes to order 2 .We have computed them to all orders in using the same techniques as those used for the two-loop amplitudes.This is by now an easy calculation so we will not discuss it further, and simply include the results in ancillary files.

Two-loop Numerical Unitarity
Our approach to the calculation of two-loop amplitudes is built on their numerical evaluation within the framework of two-loop numerical unitarity.We target independently each of the helicity amplitudes A (2,j) h with j = {0, N f }, h = {+ + +, − + +}, see eq. (2.14).The approach starts from a parametrization of the integrand of the amplitude A (2,j) h ( l ) in terms of master integrands and surface terms [25] ( l denotes the loop momenta).The master integrands are associated with master integrals, and the surface terms integrate to zero.This decomposition is naturally organized in terms of propagator structures.More precisely, we write where ∆ is the set of propagator structures Γ, P Γ is the multiset of inverse propagators ρ j in Γ, and M Γ and S Γ denote the sets of master integrands and surface terms.The coefficients c Γ,i are determined using the factorization properties of the integrand A (2,j) ( l ) in specific configurations Γ l of the loop momenta where the inverse propagators ρ j ∈ P Γ are on-shell, that is ρ j ( Γ l ) = 0 iff j ∈ P Γ .In this limit, the leading contribution to eq. (3.1) factorizes as The sum on the right-hand side is over the propagator structures Γ such that P Γ ⊆ P Γ .
On the left-hand side, T Γ denotes the set of tree amplitudes associated with the vertices in the diagram corresponding to Γ and the sum is over the scheme-dependent physical states of each internal line of Γ.We refer to eq. (3.2) as cut equations, and to its left hand side as cuts.The coefficients c Γ,i are determined by sampling the cut equations over enough values of Γ l .To construct the color-stripped products of tree amplitudes on the left-hand side of eq.(3.2), we use the unitarity-based color decomposition approach of [55,56] to include colorless particles.The dimensional-regulator dependence of cuts is determined with the approach of decomposition by particle content [49,57,58], based on dimensional reduction.We evaluate the color-stripped and dimensional-regulator-free tree amplitudes through Berends-Giele recursion [59].The resulting system of equations is then solved numerically on a given phase-space point.All numerical operations are done using finite-field arithmetic.This allows us to obtain exact coefficients c Γ,i for rational phase-space points.The latter can be generated, for instance, by using momentum-twistor variables [60].Once the c Γ,i have been determined, we obtain the decomposition of the amplitude in terms of master integrals, where m Γ,i is the master integral associated with the numerator m Γ,i ( l ).
Compared to five-parton amplitudes, we note that the fact that photons cannot be ordered leads to a proliferation of topologies.To be more explicit, consider as an example the two-loop leading-color color-ordered five-gluon amplitudes.Due to the color ordering, there are only 5 different pentagon-box topologies2 to consider, which can be labeled e.g. by the vertex attached to the gluon of momentum p 1 .For the photon amplitudes we are considering in this paper, we should multiply this number by 3!, corresponding to the possible orderings of the photons.We thus have to consider 30 different pentagon-box topologies.

Pentagon Functions
A basis of master integrals for planar five-point massless amplitudes is known [18,19].Order by order in , they can be expressed in terms of multiple polylogarithms.MPLs form a special class of functions with logarithmic singularities, and this class of functions can be equipped with algebraic structures which allow one to find relations between them [61][62][63].As a consequence, we can define a set of functions, called pentagon functions, which form a basis for the MPLs that appear in the master integrals contributing to planar two-loop five-point massless amplitudes [22,23].
After evaluating the coefficients c Γ,i in eq.(3.1), we directly obtain the decomposition of the helicity amplitudes in terms of master integrals, see eq. (3.3).Assuming that the decomposition of the master integrals in terms of pentagon functions is known, we in turn obtain a decomposition of the amplitude in terms of pentagon functions.If we denote the pentagon functions by {h i } i∈B , with B the associated set of labels, we can then decompose the amplitude as where we make explicit that two-loop amplitudes have at most poles of order −4 .The decomposition of eq.(3.4) presents a major advantage compared to the decomposition of an amplitude in terms of master integrals: it allows us to write one-loop and twoloop amplitudes in terms of a common basis of functions.It then follows from eqs. (2.16) and (2.19) that remainders themselves can be decomposed in terms of pentagon functions, 3that is, R In this paper, we will adopt the pentagon functions defined in ref. [23].Aside form the fact that they can be efficiently evaluated across phase-space, they are also defined to be a basis for the whole orbit under the symmetry group of five-point kinematics in the {1, 2}-channel (see eq. (2.4)).This is crucial for our calculation.As already mentioned previously, the amplitudes receive contributions from all orderings of the photons in the final state.An important consequence of this observation is that there is no Euclidean region for the amplitude, that is no region where the amplitude is real.

Functional Reconstruction of Master Integral Coefficients
In refs.[47,49], we argued that having two-loop corrections in mind one should reconstruct the analytic form of the coefficients r i of the decomposition of the remainders in eq.(3.5).Indeed, these are expected to be much simpler than the master integral coefficients in the decomposition of the amplitude in eq.(3.3).To better understand the simplicity of the remainder in comparison to the two-loop amplitudes for five-point QCD amplitudes, in this paper we reconstruct the coefficients c Γ,i in eq.(3.3).We thus obtain a decomposition for the amplitude that is valid to all orders in .This is the first time that such a decomposition has been made available for a two-loop five-point QCD amplitude.While this expression contains more information than required for the computation of NNLO corrections to the process in eq.(2.1), it can be used for defining the three-loop remainder required for N 3 LO corrections.For NNLO applications, we can use it to extract two-loop remainders in the form given in eq.(3.5).
Let us briefly discuss how the master-integral coefficients are computed.We proceed following the same steps as in ref. [49], which we adapt to the reconstruction of masterintegral coefficients.For simplicity of the expressions, we use a single label to index the master integrals and rewrite eq.(3.3) as where s denotes the set of five independent s ij defined in eq.(2.2) and the Levi-Civita contraction tr 5 is defined in eq.(2.3).We discuss the details of our choice of masterintegral basis in appendix A. We then make the ansatz that c i ( , s, tr 5 ) = 1 P i ( ) where we used the fact that there are no poles in that are kinematic dependent, and κ i is the maximal power of in the numerator.The polynomials P i ( ) are trivial to determine by sampling the coefficients at enough values of (for details see e.g.[36]).Since tr 5 can be written as the square root of a polynomial in the Mandelstam variables, it follows that the most generic coefficient c i,k can be written as where c + i,k ( s) and c − i,k ( s) are rational functions of the s ij .We note that which means that we can access individually each of these rational functions of the s ij by evaluating the amplitudes at parity-conjugate phase-space points.
To determine the analytic form of the rational functions c ± i,k , we numerically evaluate the master integral coefficients using the two-loop numerical unitarity method outlined in section 3.1.We recall that these numerical evaluations are performed using finite-field arithmetic which allows us to obtain exact values for the coefficients [26].We use the variables defined in appendix C of [49] to obtain a rational parametrization of phase space.This parametrization presents the major advantage of having all but one twistor variable being Mandelstam invariants, which removes any ambiguity in mapping an expression from twistor variables to Mandelstam invariants.Finally, in all calculations we set s 12 = 1 and work with dimensionless variables.The dependence on s 12 is reintroduced by dimensional analysis.
As in ref. [49], we find that the denominator of the c ± i,k ( s) can be easily determined.Indeed, we find that where the W ( s) are a subset of the so-called letters of the symbol alphabet associated with the contributing master integrals, namely the subset that is polynomial in s and not tr 5 . 4s the photons are not ordered, it is not sufficient to consider the planar alphabet which is only closed under cyclic permutations.The full non-planar alphabet [64], however, is Table 2: Characterizing data for analytic expressions of master integrals, see eqs.(3.6) and (3.7).'Max degree' is the highest polynomial degree across the different numerators c ± i,k .The column '# independent c ± i,k ' gives the dimension of the space of rational functions required to write all the c ± i,k .
closed under all photon permutations and sufficient to determine all the denominators.The powers q i,k,j are determined by evaluating the amplitudes on a one-dimensional line in phase-space [47,49].
The determination of the analytic form of the coefficients c i,k is then reduced to the determination of the polynomials n ± i,k ( s), which depend on four variables (once we set s 12 = 1).With our in-house implementation of the multivariate Newton method [36], we obtain their analytic form with coefficients in a finite field.Following the simplification procedure of the expressions described in ref. [49], we were able to reconstruct all rational coefficients from a small number of numerical evaluations in additional finite fields.We provide these results in ancillary files.
In table 2 we compile information that characterizes the complexity of the master integral coefficients.The column labeled 'Max degree' is a measure of the complexity of the functional reconstruction step.The column '# independent c ± i,k ' is a measure of the complexity of the final result.We find that the complexity of the reconstruction is comparable to the two-loop remainders of five-parton amplitudes reconstructed in ref. [49], but the complexity of the final result is much higher than in the five-parton remainder case (where the number of independent rational functions was at most O(120)).

Two-loop Finite Remainders
Once the coefficients c i in eq.(3.6) have been determined, and the expressions for the master integrals in terms of pentagon functions are known, we can obtain the decomposition of the two-loop remainder in terms of pentagon functions, see eq. (3.5).The coefficients r i have a form similar to the c i,k , that is where the r + i ( s) and r − i ( s) are rational functions of the s ij .As for the c ± i,k , the denominator of the r ± i is given by products of the letters in the subset of the (non-planar) symbol alphabet that are polynomial in the s ij .It is interesting to note that the denominator tr 2  5 is absent, even though it is polynomial in the Mandelstam invariants.This is not the case for the c ± i,k , that is there are poles in the coefficients c i ( , s, tr 5 ) that are absent in the r i ( s, tr 5 ), and the r i are thus expected to be numerically more stable than the c i .We also note that, owing to the complexity of the rational functions in c ± i,k , it is a nontrivial task to obtain the expressions for the r i from them.Finally, it is important to use the simplification procedure described in ref. [49] to obtain relatively compact expressions suitable for efficient numerical evaluation.Our expressions for the remainders are provided in ancillary files.
In table 3 we compile some data to characterize the complexity of the analytic expressions.Compared to table 2, we added the column 'Max weight' which gives a measure of the complexity of the contributing pentagon functions. 5The highest possible transcendental weight in a two-loop remainder is 4, and only one remainder saturates this bound.This number is also important for numerical evaluations, as the highest weights dominate the evaluation time of the pentagon functions.It is interesting to note that the maximal polynomial degree of the numerator we need to reconstruct is only slightly lower than in the case of master integral coefficients, see table 2, which means that the complexity of the reconstruction of remainders and of master integral coefficients is not substantially different.The complexity of the final result for master integrals is however much higher, since they depend on a larger number of independent rational functions.As expected, for NNLO applications it is thus much more efficient to work at the level of the remainders.

Reference Values
In order to facilitate the comparison with our results and to explicitly demonstrate the pole structure of the amplitudes, we present a numerical evaluation of the remainders and loop amplitudes on a randomly-chosen phase-space point in the physical region.Given that we are computing the Lorentz-invariant quantities defined in eq.(2.14), it is sufficient to specify the five independent Mandelstam variables together with the value of tr 5 , (3.12)  5: Reference evaluations of all independent bare two-loop amplitudes on the phasespace point of eq.(3.12).
The one-and two-loop amplitudes and remainders evaluated on this point are presented in tables 4 to 6. Expressions for the amplitudes in terms of master integrals for the one-and two-loop amplitudes are presented in a series of ancillary files in the directories anc/oneLoopAmplitudes/ and anc/twoLoopAmplitudes/ respectively.Expressions for the the finite remainders in terms of pentagon functions are presented in a series of files in the directory anc/remainders/.We show how to assemble these files into the full amplitudes and remainders in anc/example assembly.m,where, using an included numerical evaluation of the master integrals and pentagon functions, the expressions are combined to compute the amplitudes and remainders at the reference phase-space point (3.12) and reproduce the numbers in tables 4 to 6.

Validation
In this section we discuss the checks we performed on our setup to compute the two-loop amplitudes for three-photon production at hadron colliders.
We first computed analytic four-and five-point one-loop amplitudes.The four-point amplitudes were checked against the results of ref. [10], and we found full agreement with the quoted remainders.For the five-point amplitudes, we reproduced the results obtained from OpenLoops [65] to order 0 and numerically verified all the relations in table 1.Finally, we checked that five-point one-loop amplitudes have the correct collinear behavior up to  6: Reference evaluations of all independent one-loop and two-loop remainders on the phase-space point of eq.(3.12).order 2 , i.e., that in these limits they are either regular or factorize into products of splitting amplitudes and four-point amplitudes.
At the two-loop level, we first recomputed the four-point amplitudes, reproducing the two-loop remainders given in ref. [10].For the five-point amplitudes, we verified that they have the correct pole structure and numerically checked that they satisfy the relations in table 1.To further check the results at order 0 , we also verified that the amplitudes have the correct collinear behavior as described in appendix B.

Numerical Evaluation
Squared Finite Remainders The contribution of two-loop helicity amplitudes to physical cross sections is constructed from the finite remainders (see e.g.[52]) -specifically the squared finite remainders, summed over helicity and color states.This object, which we denote H, admits a perturbative expansion in powers of the renormalized coupling α s , which we normalize such that H (0) = 1.The O(α s ) contribution H (1) is given by and the O(α 2 s ) contribution H (2) by In this expression, the first line corresponds to the one-loop squared contribution, and the second line to the interference of two-loop and tree-level amplitudes.In the one-loop squared contributions, we have expanded the C 2 F factor for consistency with the limit in which the two-loop contributions were computed.The normalization factor in eqs.(4.2) and (4.3) is consistent with setting H (0) = 1, and in our normalization The helicity sums in eqs.(4.2) and (4.3) are performed using the relations from table 1.We note that the remainders R (2,0) +++ and R (2,N f ) +++ do not contribute to H (2) , in the same way that R (1) +++ does not contribute to H (1) .Nevertheless, R +++ does contribute to H (2)  through the one-loop squared contributions.
Numerical Evaluation Having phenomenological applications in mind, we have implemented in a C++ library the numerical evaluation of the finite remainders R )), and of H (1) and H (2) as defined in eqs.(4.2) and (4.3).This library can be obtained from a git repository [51].The library relies on PentagonFunctions++ [23] for the numerical evaluation of pentagon functions.For installation and usage instructions we refer to the README.mdfile which can be found in the root directory of the repository.
We recall that the rational coefficients in the decomposition of the remainders in terms of pentagon functions are simplified using the multivariate partial fractioning procedure outlined in [49].Furthermore, we optimize the evaluation of rational coefficients using FORM [66,67].As a result, the time spent on their numerical evaluation is negligible compared to the time spent on the evaluation of transcendental functions.On average, the evaluation time of H (2) is on the order of a few seconds per phase-space point in double precision, using the default settings of PentagonFunctions++.
Besides the evaluation speed, the calculation of the H (i) must be numerically stable.To demonstrate the stability of our results, we compare the numerical evaluation of H (2)  in double precision, which we denote H (2) double , with the evaluation in quadruple precision, which we denote H (2) quad , on a sample of 90000 phase-space points from a distribution employed in [16].We set N c = 3, N f = 5, and the renormalization scale µ R to the invariant mass of the three-photon system m γγγ .Assuming H (2) quad to be correct at least up to a relative error of ∼ 10 −16 , we define as a measure of the number of correct decimal digits in H double .In fig. 2 we show a histogram of this quantity for the 90000 sampled points on a logarithmic scale.We observe that less than 0.1% of the sampled points have an accuracy of less than four digits.This level of accuracy is more than adequate for phenomenological applications.Indeed, our implementation was already used for a Monte Carlo phase-space integration in [16], which converged to an overall integration error below 1% in NNLO differential distributions.
It is interesting to note that a good understanding of the physical properties that govern the analytic structure of scattering amplitudes can be used to explain and improve the The logarithmic distribution of decimal digits (as defined in eq.(4.5)) for 90000 double-precision evaluations of the H (2) function.The phase-space points are sampled from a distribution representative of typical phenomenological studies.
behavior of numerical algorithms.In particular, we recall that our analytic representation of the remainders (and thus of the H (i) ) only has poles that are associated with a subset of the letters of the two-loop five-point massless alphabet, see the discussion in section 3.4.While the letters corresponding to unphysical or spurious singularities can vanish inside the physical phase space, the amplitude must stay regular.This implies that large cancellations can potentially occur when the amplitude is evaluated on phase-space points that are close to the surfaces where those letters vanish.Conversely, away from these small neighborhoods, the numerical evaluations should be accurate.We observe precisely this behavior.We verified that all of the unstable phase-space points in fig. 2 are close to surfaces associated with spurious singularities.This observation provides further backing to the claim that it is important to organize the analytic structure of the rational coefficients in the expressions for amplitudes or remainders using physical considerations.For instance, once such a form is found, we can setup a robust precision-rescue system based on the fact that all the problematic regions of phase-space are explicitly known.We leave this for future work.

Conclusion
In this paper we have computed the two-loop planar corrections to the production of three photons at hadron colliders.This was achieved within the framework of two-loop numerical unitarity, coupled with analytic reconstruction techniques.Our results include expressions for both one-and two-loop amplitudes valid to all orders in the dimensional regulator.
In both cases, the amplitudes are written in terms of a set of master integrals.To our knowledge, this is the first time that master-integral coefficients have been obtained in analytic form for physical two-loop five-point scattering amplitudes.All our results are presented in a set of ancillary files.
As is well known, for NNLO (two-loop) phenomenology only the finite remainders of two-loop amplitudes are needed.By writing one-and two-loop amplitudes in terms of a basis of pentagon functions up to the required order in , we obtained a decomposition of the remainder in terms of these functions.After simplifying the coefficients in this decomposition using multivariate partial-fractioning techniques, we obtained compact expressions for the two-loop remainders of the two independent helicity amplitudes.We demonstrated that the remainders are simpler functions than the all-order amplitudes.
While over the last years there has been a number of new results for massless five-point amplitudes, when having in mind phenomenological applications is also important that the expressions for the remainders are numerically stable and can be efficiently evaluated across the relevant physical phase space.This point has been a major obstacle to computing the NNLO corrections for those processes.The expressions we obtain in this paper are ready to be used for phenomenological studies.This was demonstrated by verifying the numerical stability of the remainders when combined in a code that computes the colorand helicity-summed squared remainders.We have made this code public in a format that can be interfaced with real-radiation programs and employed to compute complete NNLO theoretical predictions.This is the first time that a two-loop massless five-point process is available in this form.Our results have already been used in ref. [16], and we expect that they will be instrumental in any future studies.index i, and the relevant scattering kinematics in terms of five-point Mandelstam invariants.This information appears in the format MI[topologyName, i, mandelstam1, mandelstam2, ...].
In this representation, the permutations act on the Mandelstam arguments and hence permutations of master integrals are codified as different Mandelstam invariants appearing in the final arguments.
Each m Γ,i denotes a numerator.These are polynomial in the loop momenta, and rational in external kinematic variables and in the dimensional regulator D. The notation used in specifying the numerators in the ancillary files is defined in table 7. There, we make use of the scalar product of the components of the loop momenta beyond 4-dimensions We collect the definitions of the master integrals employed in the one-and two-loop amplitudes in the ancillary files anc/oneLoopAmplitudes/MasterIntegralDefinitions.m and anc/twoLoopAmplitudes/MasterIntegralDefinitions.m, in the format where num is the expression of the numerator.All master integrals are "unitarity compatible", in that they have unit propagator powers, and we specify the numerator m Γ,i ( ) for a single representative of each topology.We present the external-momenta conventions and loop-momenta routing for a representative permutation of the one-loop integrals in table 8, and of the two-loop integrals in tables 9 to 11.In all cases, we take all external momenta as outgoing and we omit the loop momentum labels when they are not required to disambiguate the numerators.We can then expand both sides of (B.3), and we are particularly interested in the order (α 0 s ) 2 contributions: Ā(2) (i h i , j h j , . ..) P h (i h i , j h j ) Ā(0) (P −h , . ..) + Split (1) P h (i h i , j h j ) Ā(1) (P −h , . ..)

(B.5)
It can be useful to understand this equation graphically, as presented in figure 3.As in eq.(B.3), if the right-hand side vanishes, then the left-hand side is regular in the i||j limit.
In the context of the three-photon production amplitudes at hand, many of the collinear limits are equivalent due to the Bose symmetry of the final state.As the process is unordered, the limits can be categorized according to the particles which become collineareither (γ, γ), (q, q), (q, γ) or (q, γ).These can then be further classified according to the helicity of the involved particles.Our approach to check that our results satisfy eq.(B.5) is to evaluate both sides of the equation numerically.Note that this also tests the pentagon functions in near-singular regions of phase-space.To check their numerical behaviour, we constructed differential equations for the pentagon functions from their analytic representation, which we then solved numerically to high precision using [69].In the following, we discuss each of the different collinear configurations.Before doing so, we note that we have performed the same checks on the one-loop five-point amplitudes up to order 2 .This allows one to check the implementation of the one-loop splitting amplitudes up to O( 2 ) (see appendix B.1 below), and is also a non-trivial check of the O( ) and O( 2) contributions to the one-loop amplitudes which we have computed.(q, q) limit: This limit is not accessible in the physical region for three-photon production at hadron colliders, see eq. (2.4).For this reason equation eq.(B.5) cannot be checked in the (q, q) limit.(γ, γ) limit: The only splitting amplitude that can contribute is γ h 1 → γ h 2 γ h 3 which vanishes to all orders for any helicity.It then follows that the amplitudes should be regular when any of the photons become collinear, independently of the helicity.We have verified that this is indeed the case for our amplitudes.(q, γ) limit: For the amplitude Ā+++ , given the fact that all the photons have the same helicity, there is a single limit to check corresponding to the splitting q → γ + q.We note that

4 Figure 2 :
Figure 2:The logarithmic distribution of decimal digits (as defined in eq.(4.5)) for 90000 double-precision evaluations of the H(2) function.The phase-space points are sampled from a distribution representative of typical phenomenological studies.

Figure 3 :
Figure 3: Graphical representation of equation (B.5).White blobs represent the splitting amplitudes and hatched blobs represent the four-point amplitudes.The number of contained circles represents the loop order of the object.The types of particles are omitted from this representation, the external and internal lines representing only momenta.

Table 1 :
Relation between all helicity configurations and the two independent basis elements in eq.(2.11).P denotes a parity transformation.

Table 3 :
Characterizing data for analytic expressions of two-loop remainders, see eq.(2.19)and the discussion below for the notation.We use the same labels as in table 2. 'Max weight' is the highest transcendental weight of the pentagon functions appearing in each remainder.

Table 4 :
Reference evaluations of all independent bare one-loop amplitudes on the phasespace point of eq.(3.12).Table

Table 7 :
Notation used in specifying the numerators of master integrals.