Two-loop leading-color helicity amplitudes for three-photon production at the LHC

We calculate all planar contributions to the two-loop massless helicity amplitudes for the process $q\bar q\to \gamma\gamma\gamma$. The results are presented in fully analytic form in terms of the functional basis proposed recently by Chicherin and Sotnikov. With this publication we provide the two-loop contributions already used by us in the NNLO QCD calculation of the LHC process $pp\to \gamma\gamma\gamma$ [Chawdhry et al. (2019)]. Our results agree with a recent calculation of the same amplitude [Abreu et al. (2020)] which was performed using different techniques. We combine several modern computational techniques, notably, analytic solutions for the IBP identities, finite-field reconstruction techniques as well as the recent approach [Chen (2019)] for efficiently projecting helicity amplitudes. Our framework appears well-suited for the calculation of two-loop multileg amplitudes for which complete sets of master integrals exist.

In this work we focus our attention on QCD which is the theory most relevant to highprecision physics at the Large Hadron Collider. Within QCD, the current frontier for 2 → 1 and 2 → 2 processes is three loops [62][63][64][65][66][67][68], while for 2 → 3 scattering processes with massless partons it is two loops. A lot of work has already been carried out in this direction, mostly for planar amplitudes [69][70][71][72][73] but recently also for non-planar ones [49,57,59,[74][75][76][77]. All this progress has enabled the very recent first calculation of a 2 → 3 process at NNLO, namely, three-photon production at the LHC. This has been achieved by two separate groups [33,42] following very different methods. An essential ingredient for these calculations was the corresponding two-loop amplitude qq → γγγ. Both references have evaluated it in the leading-color approximation, following different computational approaches. The actual results for the amplitude have, until very recently, not been publicly available.
The goal of the present work is to complete this gap and present the explicit analytic result for the leading-color two-loop amplitude qq → γγγ used in the NNLO calculation of ref. [33]. We also compare our result with the one used in ref. [42] and published recently in ref. [78]. We find full agreement between the two results. Our subsequent discussion will be focused on methods used for the evaluation of the amplitude and we refer the interested reader to ref. [33] for a broader introduction to the problem and the subject, the definition of the leading color approximation for this process as well as the implications of this calculation.
The paper is organized as follows. In sec. 2 we detail the evaluation of the amplitude. Specifically, in sec. 2.1 we introduce our notation and define the finite remainder; in sec. 2.2 we describe the method for projecting helicities, while in sec. 2.3 we explain how the rational coefficients of the amplitude are derived. Our results are presented in sec. 3. They are available for download in electronic form with the arXiv submission of this work.

Notation and renormalization
We consider the partonic process where h i ∈ {+, −} denotes the helicity of the i'th parton, i = 1, . . . , 5. The indices c, c ′ denote quarks' color. All partons are massless and on-shell p 2 i = 0. Momentum conservation and on-shell conditions leave five independent parity-even Lorentz invariants s ij = (p i + p j ) 2 and one parity-odd tr 5 = 4iε p 1 p 2 p 3 p 4 . We choose the following set of variables to parameterize the amplitudes The UV renormalized amplitude for this process is denoted by where we factored out the (trivial) color dependence. We summarize the helicity configuration byh = {h 1 , h 2 , h 3 , h 4 , h 5 } and suppress the kinematic dependence for brevity. The amplitude can be expanded in α s The UV renormalized amplitude Mh is related to the bare amplitude computed in where Z q is the light quark wave-function renormalization constant. The bare coupling α 0 s is renormalized in the MS scheme according to (2.14) All UV renormalization constants are given in appendix A. The IR divergences of the UV renormalized amplitude can be factorized by means of the so-called Z operator: Once the Z-factor, the finite-remainder F and the amplitude M have been expanded in powers of α s /(4π), eq. (2.15) reduces to We define Z in the MS scheme. This completely specifies the finite remainder Fh. The explicit expansion for Z through two-loops in QCD is given in appendix A. The amplitude can be decomposed in color and electric-charge structures. For the twoloop finite remainder we find four non-vanishing contributions In the above equation, Q i are quarks' QED couplings (Q u,c,t = 2 3 e and Q d,s,b = − 1 3 e). The sum over q ′ in eq. (2.19) goes over all massless quarks. The difference between the QED couplings in the first and second lines of eq. (2.19) is due to the following: all diagrams contributing to the first line of eq. (2.19) have all three photons coupling to the external quark line (of flavor q) while the diagrams contributing to the second line of eq. (2.19) have one photon coupling directly to the external quark line (of flavor q) and two photons coupling to an internal fermion loop with a flavor q ′ . Example diagrams for the different contributions can be found in fig. 1. In this work we do not consider diagrams with massive quark loops. We also note that the contribution proportional to C F q ′ Q 3 q ′ from diagrams with three photons coupling to an internal fermion loop (see fig. 1) vanishes by Furry's theorem.
The leading color contribution of Fh (2) (qq → γγγ) is proportional to N 2 c and, as follows from eq. (2.19), is only dependent on a linear combination of the first two factors in the first line of eq. (2.19): The phenomenological analysis in ref. [33] is based on eq. (2.20). The justification for the use of this approximation can be found in that reference.
Despite recent progress, the non-planar diagrams in this process are still beyond reach. For this reason, in this work we derive and present in analytic form the planar results for the following two factors:

Helicity projections
The process (2.1) has two independent helicity amplitudes. As such we choose All other helicities can be obtained from this set by conjugation and/or permutation of external momenta. Since at tree-level F (0),h − = 0, the one-loop contribution to this helicity is finite and its two-loop correction does not contribute to the squared matrix element for this process through two loops. In order to extract the helicity amplitudes (2.22) in the 't Hooft-Veltman scheme we employ the projection method proposed in ref. [79]. Similar approaches have been advocated in refs. [80,81]. The essence of this method is that, for a given helicity, it provides explicit prescription for constructing the external wave-functions. For the process under consideration, these wave-functions can be factored out by writing the amplitude in the following way The construction of the external wave-functions is done in 4 dimensions. On one hand this significantly simplifies the construction of a basis of polarization vectors. On the other, it introduces scheme dependence into the bare loop amplitude. As demonstrated in refs. [79,82,83], however, this scheme dependence does not affect the finite remainder of the amplitude in the limit ε → 0.
Our construction of the fermionic wave functions introduces the matrix γ 5 . The fact that wave functions are constructed in 4 dimensions implies that we only use 4-dimensional identities for γ 5 . With the help of these identities we eliminate all occurrences of γ 5 and, in this process, trade them for objects involving ε µνρσ . Multiplying (and dividing) the term proportional to ε µνρσ by (tr 5 ) 2 and recalling that (tr 5 ) 2 is a polynomial function of the invariants s ij (see eq. (2.9)), we eliminate all occurrences of ε µνρσ using and then promoting all indices to d dimensions. Once this has been achieved, the subsequent contractions of external polarization states with the rest of the amplitude are performed in d dimensions. The remaining factor ε p 1 p 2 p 3 p 4 is converted into tr 5 which is then treated as an independent kinematic variable that is included in the rational coefficients A b defined in sec. 2.3 below. We next derive explicit expressions for the wave functions of the final state vector particles of given helicity As a first step we replace those with new linear polarization vectors polarized along two directions X and Y . The polarization vector along X is defined through the following ansatz As a reference vector for all vectors ε µ i,X we chose the vector q µ = p µ 1 +p µ 2 . The coefficients c X i,n , n = 1, 2, 3 are determined from the system of normalization and orthogonality conditions for the polarization vectors: (2.28) The polarization vector along the direction Y is given by where we have used the conventions of ref. [79] for the Levi-Civita symbol: ε 0123 = +1 and ε 0123 = −ε 0123 . The normalization factors N i,Y are determined from the condition ε 2 i,Y = −1. Lastly, we note that the above construction does not fix the vectors ε i,X and ε i,Y uniquely. There are two possible solutions which correspond to the change ε i,X/Y → −ε i,X/Y . The overall signs of these two vectors are chosen in such a way that the vectors p i , ε i,X and ε i,Y form a right-handed coordinate system.
The fermion wave functions are treated in the following way. The spinor part of the amplitude (2.23) has the following structure The matrix u ⊗v in the above equation can be rewritten in the following way: for some matrix N , to be specified below, and N ≡ūN v = 0. The outer spinor products read The matrix N depends on the process-specific kinematics. In particular, it is linearly independent of p 1 and p 2 , otherwiseūN v = 0. Since for helicity configurations with h 1 = h 2 the amplitude (2.1) vanishes to all orders, in the following we only consider the case h = h 1 = −h 2 . For this helicity configuration the matrix N is given by (2.34) For this choice of N , eq. (2.31) takes the form

Reduction to pentagon functions
In this work we consider four independent amplitude structures. They correspond to the two color factors in eq. (2.21) and the two helicities in eq. (2.22). We would like to express these structures in terms of transcendental functions and transcendental constants with rational coefficients. Since we assume that the set of functions forms a basis, the evaluation of an amplitude structure is equivalent to the evaluation of all its rational coefficients. To this end, we have built an automated framework that uses finite-field methods to numerically evaluate, interpolate and then reconstruct the exact analytical expressions for these coefficients. Throughout this section, the term numerical should be understood to refer to finite-field numerics. We will now provide a detailed account of our framework.
Any bare scalar 2-loop amplitude, or amplitude structure, M is expressed as a linear combination of 2-loop scalar integrals I b : (2.36) The coefficients A b are rational functions of the kinematic invariants s ij , defined in sec. 2.1, and polynomials in ε. As explained in sec. 2.2 the coefficients A b are also linear functions of the parity-odd kinematic variable tr 5 . The planar two-loop integrals I b appearing in the amplitude of interest in this paper are defined [43] as follows: where the propagators Π i are defined in table 1 and n i ≡ n i (b) ∈ Z. These 11 propagators form a complete basis of bilinears through which any scalar numerator structure in an integrand can be expressed. A maximum of 8 propagators appear in the denominator of any given integrand; for the three remaining propagators, denoted spurious propagators, the corresponding indices n i satisfy n i ≤ 0. We can accordingly classify the integrals into two topologies, which are represented pictorially in fig. 2. In the C 1 topology, propagators 7, 10, and 11 are spurious. In the C 2 topology, propagators 6, 7, and 10 are spurious. Following the standard Integration-By-Parts (IBP) approach, we map the integrals I b onto a basis of master integrals, M c (2.38) The coefficients B bc are rational functions of s ij and ε. All coefficients required in this calculation are known analytically from ref. [43]. The IBP identities map all required C 1 integrals onto a basis of 61 master integrals 1 , and all required C 2 integrals onto a basis of 28 master integrals. The 28 C 2 master integrals can themselves be identified with C 1 integrals, in some cases by permuting the external momenta {p j }. The next step in our calculation is to map the master integrals M c onto a basis of transcendental functions and constants. In this work we choose the basis of functions provided recently in ref. [59]. That reference solves a set of specially designed loop master integrals where the basis {t e } is built out of sums and products of the transcendental functions F and transcendental constants tcr and tci that are defined in ref. [59] as well as the unity. The coefficients D de are rational functions of the kinematic invariants s ij and linear functions of tr 5 . They also depend on ε, however, they are only known as a series expansion of sufficient depth appearing in eq. (2.47) are defined as The coefficients B can be split into two parts: one which is proportional to tr 5 and one which is independent of it. Both parts are rational functions of the invariants s ij . The reason only the first power of the parity-odd variable tr 5 appears in the final result is that (tr 5 ) 2 is itself a polynomial of s ij , see eq. (2.9).
The calculation of the coefficients G (n) e is based on eq. (2.47) and proceeds as follows. We first note that the coefficients A e . The problem with this strategy is that the size of the rational expressions that need to be combined becomes huge which hampers their subsequent simplification. Such a strategy was followed by us in ref. [33] for the evaluation of the squared amplitude for qq → γγγ. We refer the reader to that reference for more details about the subtleties of such an approach.
In this work we follow an alternative strategy for the evaluation of the coefficients G e are then passed to the FireFly library [85] in order to reconstruct the exact analytical expressions for the coefficients G (n) e . We have created an automated framework which is designed to calculate simultaneously multiple amplitudes (that have the same kinematics) following the finite-field evaluation approach just described. The reason it may be advantageous to compute several amplitudes at the same time is efficiency, noting that only the coefficients A To fully specify our approach we need to describe one more feature which has to do with momentum crossings. The coefficients A do not require any further crossing since all crossings are already included at the diagrammatic level. Similarly, the functions {t e } are already defined in such a way that all possible crossings have been already implemented in their definition [59]. This is one significant difference with respect to the functional basis constructed in ref. [47] and used by us in ref. [33]. The basis of ref. [47] does not include momentum crossings and if they are required the user needs to implement those. As a result of such crossings one generally arrives at a functional basis which is non-minimal. To reduce the extended set of functions to a minimal set, functional identities need to be derived and applied; see ref. [33] for more details on this point.
Where momentum crossings still need to be applied is the part of eq. (2.47) that involves IBPs. Specifically, this affects the coefficients B  fig. 2. In principle any crossing of the IBP solution can be obtained analytically by rearranging the solution in terms of the crossed kinematic invariants. Such a strategy would be impractical in our numerical approach due to the significant size of the IBP solutions and the very large number of crossings.
The way we deal with crossings in our practical implementation is as follows. The master integral reductions E for all amplitudes M and we only store running 2 We use a multi-core computing cluster. Each computing thread is assigned independent finite-field points and evaluates all the amplitudes at the assigned points, one point at at time.

Results
Following the approach described in the previous section, in this work we calculate in analytical form the two helicities eq. (2.22) of the finite remainder for the process qq → γγγ in the 't Hooft-Veltman scheme. In this work we have not included any non-planar contributions i.e. we work in the approximation eq. (2.20). Factoring out the phase-dependent part of the leading order amplitude we write the reconstructed finite remainder for the helicityh + defined in eq. (2.22) as: For the helicityh − the tree-level amplitude vanishes and we write The functions Rh (i) , i = (1, N 2 c , C F n f ), have the following structure We have checked that the above result agrees numerically with our previous calculation performed in ref. [33]. This comparison does not include the terms ∼ n f since those were not computed in ref. [33]. Given the two calculations were performed with almost completely independent methods this represents a strong check on the correctness of eqs. (3.1,3.2). Since one of the main objectives of the present work is to document the calculation of the two-loop amplitude used in ref. [33], we will next explain in some detail how the two calculational approaches differ from each other.
Ref. [33] computed directly the squared amplitude, while in this work we compute the helicity amplitudes and the squared amplitude is obtained by squaring and crossing them numerically. Ref. [33] used the polygon functional basis of ref. [47] while in the present work we use the basis of ref. [59]. The main difference between the two functional bases was explained in the previous section. Our calculation in ref. [33] was fully analytic while the current calculation uses finite-field numeric evaluation and reconstruction techniques to obtain the rational coefficients. Besides the differences already mentioned, the two results differ significantly as far as evaluation times are concerned. The main reason for this is the difference in evaluation times and numerical precision between the two functional bases. The larger size of the result in ref. [33] is immaterial given the time needed for the evaluation of the most complicated functions in that basis.
We have also checked that our results eqs. (3.1,3.2) agree with the results in ref. [78]. Since we have computed a different helicity combination relative to the helicities published in ref. [78], a direct analytic comparison between the two is complicated by the fact crossing of external legs is required. We have analytically checked most structures which are simple enough to cross, while the complete expressions (for each helicity and color factor) have been compared numerically with high numerical precision (30 digits) and full agreement between the two calculations has been found in all cases. There are certain differences in the way the present calculation and the one in ref. [78] were performed that make the agreement between the two calculations highly nontrivial. First, the generation of diagrams is based on different approaches. Furthermore, in this work we use the analytic solutions of the IBP equations from ref. [43] while ref. [78] solves the IBP equations numerically for each point in which the amplitude is evaluated and being reconstructed. Most importantly, in this work we use an alternative approach to the projection of helicity amplitudes which is very different from the one employed in ref. [78].

Conclusion
In this work we calculate the planar contributions to the two-loop helicity amplitude for the process qq → γγγ. The result is presented in fully analytic form and is available for download in electronic form with the arXiv submission of this paper. This result, written in an alternative form, was used in the first NNLO calculation of a 2 → 3 LHC process [33].
The helicity amplitudes are expressed in terms of the functional basis of ref. [59] which allows fast and efficient numerical evaluation of the amplitudes. These functions' speed of evaluation is sufficiently high to allow the direct use of this amplitude in the calculation of NNLO cross-sections without the need for intermediate interpolation.
In our calculation we have utilized Chen's recently proposed approach [79] for the efficient projection of helicity amplitudes in multileg/multiloop processes. We have found the approach very easy to implement and use especially since it does not involve investigations of tensor bases that grow with the number of loops.
We have found complete agreement between our calculation and the recent independent calculation of the same amplitude in ref. [78]. Such an agreement represents a highly nontrivial check on both calculations given they use very different computational approaches.
Our calculation is derived from an automated framework we have created for the calculation of generic two-loop massless 5-point gauge theory amplitudes. We hope it will prove suitable for many other potential applications.

Acknowledgments
We acknowledge helpful discussions with Long Chen about the helicity projection method of ref. [79] and with Fabian Lange and Jonas Klappert on the use of the library FireFly

A Renormalization constants
The UV renormalisation constant of the quark wave function through order O(α 2 s ) reads while the renormalization constant Z αs = Z 2 g up to power α s (higher powers are not required since the tree-level amplitude for pp → γγγ has a zero power of α s ) is given by Z g = 1 + 1 ε α s 4π 4(n l + 1)T F − 11C A 6 . (A. 2) The contribution from heavy flavours is not considered in this work, therefore we have Z q = 1 and the term (n l + 1) term in Z g reads n l . For the IR renormalization of the amplitude for the process qq → γγγ the color-space matrix Z is needed. Through order O(α 2 s ) it is given by with the abbreviation l µ = log −µ 2 /s 12 .