Three-loop helicity amplitudes for diphoton production in gluon fusion

We present a calculation of the helicity amplitudes for the process $gg\to\gamma\gamma$ in three-loop massless QCD. We employ a recently proposed method to calculate scattering amplitudes in the 't Hooft-Veltman scheme that reduces the amount of spurious non-physical information needed at intermediate stages of the computation. Our analytic results for the three-loop helicity amplitudes are remarkably compact, and can be efficiently evaluated numerically. This calculation provides the last missing building block for the computation of NNLO QCD corrections to diphoton production in gluon fusion.

from NNLO, the gg partonic channel opens up. There are two contributions to this: tree-level corrections of the form gg → γγ + qq and loop-induced corrections gg → γγ. Phenomenologically, the former are typically very small and we will not discuss them further. The loop-induced contribution is instead quite interesting. First, the large gluon flux at the Large Hadron Collider (LHC) compensates for the α s suppression, making it important for precision phenomenological studies. Being a new channel, it has all the features of a leading order process, in particular large perturbative uncertainties. Moreover, being gluon induced one expects particularly large radiative corrections. This has spurred many investigations, which upgraded the precision in this channel to next-to-leading order (NLO) in QCD, i.e. to O(α 3 s ) [18,19]. Given the ever-increasing experimental precision on diphoton measurements [20,21], it becomes interesting to try and push the theoretical precision even further and consider NNLO corrections to the gg → γγ process. While this is desirable for a variety of LHC analyses, it is of particular importance for Higgs studies. Indeed, in this case there is a subtle signal/background interference effect between the gg → H → γγ signal and the continuum gg → γγ background, which is known to modify the Higgs line-shape [22]. This effect can in turn be used to constrain the Higgs boson total decay width [23]. This kind of investigations require an exquisite experimental control, see e.g. ref. [24], as well as robust control of theoretical predictions for both the signal and the background processes. Several in-depth analysis [25][26][27][28] suggest that reaching NNLO QCD accuracy in the gluon channel is desirable. A major step towards the calculation of full NNLO QCD corrections to gg → γγ has been made very recently with the computation of NLO QCD corrections to the gg → γγ + j process [10,11]. In this paper, we present a calculation of the last missing ingredient, the three loop virtual amplitude for the gg → γγ process.
The rest of this paper is organised as follows. In sec. 2 we set up our notation and discuss the generic kinematics features of the gg → γγ process. In sec. 3 we briefly review the approach of refs [29,30] to the calculation of helicity amplitudes that we adopt here. In sec. 4 we provide more technical details on our three-loop calculation. In sec. 5 we discuss the ultraviolet and infrared structure of the scattering amplitude, and define the renormalised finite remainders which are the main result of this paper. In sec. 6 we document the checks that we have performed on our calculation, and briefly describe the general structure of our result. We also present analytic formulas for the three loop finite remainder for the simplest helicity configuration. The analytic formulas for all the relevant helicity configurations can be found in computer-readable format in the ancillary material that accompany this submission. Finally, we conclude in sec. 7.

Notation and kinematics
We consider virtual QCD corrections to the production of two photons through gluon fusion mediated by light quarks. The signs of the momenta are chosen such that all momenta are incoming, p 1 + p 2 + p 3 + p 4 = 0. All particles in the process are on the mass-shell, The kinematics is fully described by the usual Mandelstam invariants In the physical scattering region, one has s > 0, t < 0, u < 0. For later reference, we also introduce the dimensionless ratio where 0 < x < 1 in the physical region. We work in d = 4−2 dimensions to regulate ultraviolet (UV) and infrared (IR) divergences. To be precise, we adopt the 't Hooft-Veltman (tHV) scheme [31], i.e. we perform computations for generic d but we constrain all the external particles and their polarisations to live in the physical d = 4 subspace. This allows us to simplify the calculation compared to the Conventional Dimensional Regularisation (CDR) case, where internal and external degrees of freedom are treated as d-dimensional.
We write the scattering amplitude for the process in eq. (2.1) as where a j is the colour index of the gluon of momentum p j and j,µ (p j ) is the polarisation vector of the vector boson of momentum p j . For convenience, we have extracted the leading-order electroweak coupling written in terms of the fine structure constant α, where e = √ 4πα is the unit of electric charge. We are interested in the QCD perturbative expansion of eq. (2.4) where α s = α s (µ) is the MS renormalized QCD coupling and the superscript indicates the number of loops L. We find it convenient to express the result for A (L) in terms of the quadratic Casimir invariants of theory C A and C F . They are defined through where f abc and T a ij are the SU (3) structure constants and the generators in the fundamental representation, respectively. We normalise the generators as In QCD, C A = 3 and C F = 4/3. The Feynman diagrams for the process eq. (2.1) can be naturally separated according to whether the two photons couple to the same or to two different closed fermion lines. We then introduce the following short hands for the respective electromagnetic coupling structures where the sums run over n f light quarks and Q f is their charge in units of e, i.e. Q u,c = 2/3, Q d,s,b = −1/3. For QCD with 5 flavours, the structures in eq. (2.8) evaluate to (n V f ) 2 = (1/3) 2 = 1/9 and n V 2 f = 11/9.

The helicity amplitudes
In this section, we explain how one can efficiently calculate the amplitude in eq. (2.4) for specific helicities. We start by discussing the tensor A µνρσ . It can be expanded as where F i are scalar form factors 1 and Γ µνρσ i are independent tensor structures constructed using external momenta {p µ i } and the metric tensor g µν . With three independent external momenta, the total number of tensor structures that one can write is 138, see e.g. [32]. Since A µνρσ has to be contracted with the external polarisation vectors µ i , one can use the physical conditions p i · i = 0 to remove all tensors proportional to p µ 1 , p ν 2 , p ρ 3 , p σ 4 . This removes all but 57 structures. By making a specific choice for the reference vectors of the external gauge bosons, one may eliminate further redundancies. A convenient choice is to impose i · p i+1 = 0 , where i = 1, ..., 4 and p 5 ≡ p 1 . (

3.2)
This leaves one with 10 independent structures, that we choose as For notational convenience, we define the 10 independent structures T i = Γ µνρσ i 1,µ 2,ν 3,ρ 4,σ (3.4) and refer to them, with a slight abuse of language, as tensors. The scattering amplitude eq. (2.4) can then be written as We stress that eq. (3.5) is valid at any perturbative order and for any space-time dimension.
In four dimensions, it is easy to see that only 8 out of the 10 tensors T i are actually independent. It turns out that in the tHV scheme, it is possible to separate the purely four-dimensional tensor structures from the −2 -dimensional ones through a simple orthogonalisation procedure [29,30]. We briefly sketch how this can be done for our process, and refer the reader to refs [29,30] for a thorough discussion. Following ref. [30], we introduce a new tensor basis T i where the first 7 tensors are identical to the ones introduced before while T 8 is a symmetrised version of T 8 It turns out [30] that these 8 tensors span the physical d = 4 subspace and do not have any component in the −2 directions. The last two tensors T 9,10 can then be chosen in such a way that they are constrained to live in the −2 subspace. This can be achieved by simply removing from the original T 9,10 their projection along T 1...8 where the projectors P i are defined through The explicit form of the P i projectors relevant for our case can be found in ref. [30]. The new tensors T 9,10 read (3.11) The tensors T 9,10 so constructed identically vanish if they are computed using physical polarisation vectors in d = 4 space-time dimensions and can be safely dropped if one is after tHV helicity amplitudes [30]. For a given helicity configuration we then write where T i,λ 1 λ 2 λ 3 λ 4 are the tensors evaluated with polarisation vectors for well-defined helicity states λ i . It should not be surprising that the generic helicity amplitude can be parametrised in terms of 8 independent structures. Indeed, in four dimensions we would need to consider 2 4 = 16 independent helicity amplitudes. However, half of them can be related by parity, which leaves us with 8 independent helicity states. These are in one-to-one correspondence with the 8 form factors F i .
When dealing with helicity amplitudes, we find it convenient to factor out a spinor function carrying the relevant helicity weight. We achieve this by writing and We note that we have chosen the spinor functions in eq. (3.14) following ref. [18]. The expressions for the spinor-free amplitudes f λ 1 λ 2 λ 3 λ 4 can be easily obtained by computing the relevant T i with polarisation vectors for fixed helicity states. We also note that we define "±" helicity states as 2 where q j is the reference vector for the boson j, irrespective of whether the particles are in the initial or the final state.
We have written eqs (3.14,3.15) for only 8 helicity states. The 8 remaining ones can be obtained from these by exploiting parity invariance, where −λ i indicates the opposite helicity of λ i . We also note that the helicity amplitudes must obey Bose symmetry, i.e. they must be symmetric under the exchange of 1 ↔ 2 and/or 3 ↔ 4. In terms of the spinor-free amplitudes, this implies These relations provide non-trivial checks for our results.

Details of the calculation
The spinor-free helicity amplitudes f λ 1 λ 2 λ 3 λ 4 can be computed as perturbative series in the QCD coupling constant α s . For a generic helicity configuration we introduce the short hand λ = (λ 1 , λ 2 , λ 3 , λ 4 ) and write where α s,b is the bare strong coupling constant and f is the bare perturbative coefficient of the helicity amplitude. Since the leading order contribution f to the production of two photons in gluon fusion already involves one-loop integrals, the next-to-next-to-leading order contribution f As explained in sec. 3, we can obtain the helicity amplitudes by computing the F i , i = 1, ..., 8 form factors. In principle, this can be achieved straightforwardly by applying the projectors P i , i = 1, . . . , 8, of sec. 3 to the sum of all the relevant Feynman diagrams. At three loops, this leads to a sum of terms of the form where k i , i = 1, 2, 3, are the loop momenta, D i are the propagators of the graphs and n i are non-negative integers. Following previous work [34,35], the integration measure for every loop is defined as It is convenient to treat propagators and scalar products involving the loop momenta on the same footings. We do this by writing scalar products in the numerator as additional propagators raised to negative powers. For our problem, there are 6 scalar products of the form k i · k j and 9 of the form k i · p j , so we can write a generic Feynman integral of the form eq. (4.2) as where now n i can also be negative integers. We refer to each set of inequivalent {D 1 , ..., D 15 } as an "integral family". Within each family, it is well known that not all the integrals are linearly independent. Indeed, Feynman integrals satisfy integration-by-parts (IBP) identities [36] of the form where v j can be any loop or external momentum. In principle, it is possible to use these identities to express all the F i form factors in terms of a minimal set of independent "master integrals" (MI) [37]. While all the steps described above are well-understood in principle, the complexity involved in intermediate stages grows very quickly with the number of loops and external scales. In our case, the three-loop calculation involves 3 different families, each of which can contribute with 6 independent crossings of the external legs, and more than 4 × 10 6 integrals to the amplitude. Moreover, using (4.5) directly would lead to a very large number of equations involving also many additional auxiliary integrals. We now describe the procedure that we have adopted to keep the degree of complexity manageable. First, we generated all Feynman diagrams with Qgraf [38] and mapped each diagram to an integral family using Reduze 2 [39,40] to generate the required shifts of loop momenta. At this stage, it is useful to group diagrams that present similar structures together and perform the P 1,...,8 projections for each of these groups separately. This can be done by keeping together diagrams that can be mapped to the same crossing of the same integral families. This allows us to reduce redundancy in the algebraic manipulations required. Examples of top sectors from our three families of integrals are depicted in Fig. 1, while their complete definition can be found in the ancillary files. To evaluate the contributions to the form factors, we performed the colour, Lorentz and Dirac algebra as well as further symbolic manipulations described in the following with Form [41].
We find it important to stress that by expressing the result for each F 1,...,8 in terms of a minimal set of integrals under crossings and shift symmetries, prior to performing the actual IBP reduction, we saw a significant decrease in complexity. This is expected, as many equivalent integrals are combined together and redundant structures are removed. After this simplification, we used the F i to construct the spinor-free helicity amplitudes, and collected the contributions to different colour structures. This way, we arrived at a minimal set of gauge-independent building blocks. We found it useful to partial fraction the rational functions with respect to x in order to reduce their complexity.
The next step is the actual IBP integral reduction. We did this using an in-house implementation of the Laporta algorithm [37], Finred [42], which exploits syzygy-based techniques [43][44][45][46][47][48] and finite-field arithmetics [42,[49][50][51]. We found in this way that the three loop helicity amplitudes can be expressed in terms of the 221 MIs computed in ref. [34] and crossed versions of them, for a total of 486 MIs. 3 We stress that these MIs are pure functions, i.e. they do not have any non-trivial rational functions of x or d as prefactors. Before inserting the IBP relations into the amplitude, we partial fractioned them with respect to both d and x. We found that this step is crucial to keep the complexity under control. Finally, we performed one last partial fraction decomposition of the full amplitude after we wrote it in terms of MIs.
As a last step, we expand in and substitute the analytic results for the MIs. All of the integrals required for our calculation were computed in ref. [34]. Their expansion can be written in terms of Harmonic Polylogarithms (HPL), that we define iteratively as 4 G(0, . . . , 0 n times ; x) ≡ ln n x n! , G(a n , ..., a 1 ; x) = x 0 dz z − a n G(a n−1 , ..., a 1 ; z), (4.6) with G(x) = 1 and a i ∈ {0, 1}. For our case, we need to consider polylogarithms up to weight 6, i.e. n = 6 in eq. (4.6). We used the Mathematica package PolyLogTools [54] to manipulate HPLs up to weight 5, augmented by a straightforward generalisation of its routines up to the required weight 6, as well as an independent package for multiple polylogarithms written by one of us. As expected from the fact that there are fewer weight ≤ 6 HPLs than MIs, we observed a noticeable decrease in complexity for the amplitude after we expressed it in term of HPLs. We summarise the degree of complexity of the various steps discussed above in Tab. 1. Before presenting our results, we note that although the MIs have been computed in ref. [34], for this calculation we have decided to recompute them as an independent check. We used the same definitions for the MIs as ref. [34], and followed the same strategy outlined in that reference for obtaining their analytic. First, since the basis [34] is pure and of uniform weight [55] the MIs obey very simple differential equations where M is a vector whose components are the MIs and A i are constant matrices. Using the basis of ref. [34] we have rederived the differential equation from scratch and found 3 We note that while the reductions of some integrals were already known from earlier calculations [3,35], for this process we had to reduce a significant number of new integrals compared to those references. 4 Note that we use the GPL notation of ref. [52], rather than the original HPL notation of ref. [53]. agreement. Given the simple form of eq. (4.7), it is straightforward to iteratively solve it order by order in , modulo boundary conditions. The only non-trivial issue is how to fix the latter. Very interestingly, the authors of ref. [34] noted that at three loops it is enough to impose regularity conditions to fix all boundary conditions, apart from one simple overall normalisation. The main idea is to look at the differential equation near singular points s → 0, t → 0, u → 0. Let us consider s → 0 as an example. In this limit the general solution of eq. (4.7) behaves like M ∼ s As M 0,s , where M 0,s is a constant vector. It was argued in ref. [34] that the MIs considered here can only develop branch cuts of the form s −α with α > 0. This implies that the coefficient of s α in s As M 0,s must vanish for α > 0. As a consequence, there must exist non-trivial relations between different MIs in the s → 0 limit. When combined with analogous relations derived from the limits t, u → 0, the authors of ref. [34] found that for the case under study one can completely constrain all the boundary conditions up to an overall normalisation factor. We have independently verified that this is the case, which allowed us to rederive an analytic expression for all the master integrals. We have then verified that our results to weight 6 are identical to the ones of ref. [34], provided that the latter are analytically continued to the physical Riemann sheet. Since in ref. [34] final results are only presented for one single crossing, for convenience we decided to provide analytic results for all the three-loop master integrals and all their crossings in the ancillary files accompanying this publication. We also provide weight 6 results for a uniform-weight basis of the two-loop integrals.

UV renormalisation and IR regularisation
Following the steps outlined above, we obtained analytical expressions for the bare helicity amplitudes f are affected by both ultraviolet (UV) and infrared (IR) divergences, which manifest themselves as poles in the dimensional regularisation parameter = (4 − d)/2. While the former are removed by UV renormalisation, the latter can be regularised using universal IR operators acting on lower-loop amplitudes. We now discuss in detail how this can be done.
We first consider UV divergences. We define α s (µ) to be the renormalised strong coupling constant in the MS scheme at the scale µ The first two coefficients of the QCD beta function read We then expand the spinor-free helicity amplitudes f λ in terms of the renormalised strong coupling α s (µ) as The expression for the renormalized amplitudes f (L) λ can be obtained by substituting eq. (5.1) in eq. (4.1) and expanding in the renormalised coupling. For convenience, we will set µ 2 = s in the following. The result for arbitrary scale can be easily obtained using renormalisation group methods.
We now consider IR divergences. The IR structure of the amplitude is governed by the soft and collinear behaviour of virtual quarks and gluons and it is universal, i.e. it only depends on the colour and nature of the external legs. This allows one to write the renormalised amplitude as where f (i,fin) λ are finite in four dimensions. The IR structure is encoded in the operators I i , that for our case (with µ 2 = s) read [56] where K is the next-to-leading-order coefficient of the cusp anomalous dimension and [57] H g = 1 2 In eqs (5.5) we used the fact that diphoton production in gluon fusion starts at one loop. The finite remainders for the helicity amplitudes f (L,f in) λ up to three loops are the main result of this paper, and we provide analytic results for them in the ancillary files.

Checks and structure of the result
We have performed various checks on the correctness of our results. First, we have employed two derivations of the three-loop F i form factors at the integrand level and verified that they agree. We have also compared the one-and two-loop helicity amplitudes against the results of ref. [18] and found agreement. To validate our numerical evaluation procedure, we also checked the helicity-summed one-loop squared amplitude against OpenLoops [58,59], and one helicity configuration at two loops against MCFM [60,61]. Finally, we have verified that the UV and IR poles up to three loops follow the structure described in the previous section. This provides a strong check of the correctness of the three-loop amplitudes.
We now discuss the general structure of our result. The amplitude can be expressed in terms of the two quadratic Casimirs C A and C F and the flavour structures n f , n V f and n V 2 f defined in eq. (2.8). At L loops the amplitude is a homogeneous degree-L polynomial in these 5 variables. At one-loop, the amplitude is only proportional to n V 2 f , since the two photons must both couple to the same fermion line. At two loops, the structures n V 2 f × {C F , C A } appear in the bare amplitude. The finite remainder contains in addition a term proportional to n V 2 f n f stemming from β 0 in the UV/IR regularisation. We note that there is no (n V f ) 2 contribution. It is easy to understand why this is the case. The (n V f ) 2 colour factor only appears if the two photons are attached to two different (closed) fermion lines. Such diagrams do appear at two loops, but they are of the form of two γgg * one-loop triangles connected through a gluon propagator. Due to an argument analogous to Furry's theorem, these diagrams give no net contribution to the amplitude. A similar argument allows one to conclude that there is no net contribution from Feynman diagrams with colour factors n f (n V f ) 2 at three loops. Furthermore, it is easy to see that the structure n 2 f n V 2 f is absent in the three-loop bare amplitude. 5 Since there is no (n V f ) 2 contribution at lower loops, the (n V f ) 2 term in the bare three loop amplitude must be finite. We observe, however, that it is non-zero. Indeed, at three loops this colour factor appears in triple-box diagrams for which the Furry argument outlined above is no longer applicable.
We now move to the discussion of the kinematic features of the three-loop amplitude, i.e. its x dependence. The amplitude contains terms of the form G(a 1 , ..., a n ; x)/x k (−2 ≤ k ≤ 2) and G(a 1 , ..., a n ; x)/(1 − x) k (1 ≤ k ≤ 2) , where a i ∈ {0, 1}, 0 ≤ n ≤ 6, and G are the Harmonic Polylogarithms defined in eq. (4.6). Instead of the HPLs, we 5 We note however that the n 2 f n V 2 f structure contributes to our finite remainders, since it is induced by the n f dependence of the UV/IR counterterms.
found it useful to also consider the alternative functional basis described in ref. [35] to speed up the numerical evaluation of the final result. Using the algorithm of [62], we have constructed a basis of logarithms, classical polylogarithms and multiple polylogarithms to rewrite the HPLs without introducing any new spurious singularities. We used products of lower weight functions whenever possible and preferred functions whose series representation requires a small number of nested sums. In this way, we found that 23 independent transcendental functions and products thereof suffice to represent our HPLs up to weight 6.
In the ancillary files, we provide our analytic results written both in terms of HPLs and in terms of this minimal set of functions. For convenience, we also provide results for the finite remainders of the one-and two-loop helicity amplitudes up to weight 6. Finally, we present our results. Although intermediate expressions are rather complicated, see Table 1, we find that the final results are remarkably compact. The λ = (++++) helicity configuration is particularly simple. This is of course expected, since the one-loop amplitude does not have support on any cut, hence it is purely rational rather than a weight-2 function. This simplicity persists at higher loops. For illustration, we now report here the result for the finite reminders defined in eq. (5.5) up to three loops for this helicity configuration. At one and two loops one has f (1,fin) At three loops, we write the finite remainder as f (3,fin) iπ + 43 9 L 0 x − 7 9 x 2 (L 0 − L 1 ) 2 + π 2 , 16 3 L 0 x + 4 3 x 2 (L 0 − L 1 ) 2 + π 2 , L 0 x + 1 36 iπ . (6.5) In eq. (6.5), we have defined L 0 = ln(x) and L 1 = ln(1 − x). We note that these are the only transcendental functions that are needed to describe our result. Although the results for the remaining helicity configurations are still rather compact, they are much larger than for the λ = (+ + ++) case. We provide them in electronicallyreadable format, attached the arXiv submission of this paper. In Figure 2 we plot our result for the one-, two-and three-loop finite remainders M We fix α s = 0.118 and show graphs for the helicity configurations λ = (++++), (−+++), (+ + −+), (− − ++) and (− + −+). All the other helicity amplitudes can be obtained from these through Bose symmetry (x ↔ 1 − x) and parity.

Conclusions
In this paper, we have computed the helicity amplitudes for the process gg → γγ in threeloop massless QCD. This is the last missing ingredient required for the calculation of the NNLO QCD corrections to diphoton production in the gg channel. For our analytical threeloop calculation, we have adopted a new projector-based prescription to compute helicity amplitudes in the 't Hooft-Veltman scheme. The expressions at the intermediate stages of our calculation were quite sizable, and we employed recent ideas for the demanding integration-by-parts reductions. Our final results though are remarkably compact. They can be expressed either in terms of standard Harmonic Polylogarithms of weight up to six, or in terms of only 23 transcendental functions defined by up to three-fold sums. This makes the numerical evaluation of our result both fast and numerically stable. Analytical results for both choices of the transcendental functions are provided in the ancillary files that accompany this publication.
We envision several possible future directions of investigations. On a more phenomenological side, it would be interesting to combine our results with those of refs [10,11]  obtain NNLO predictions for the gg → γγ process. On a more theoretical side, the simplicity of our final results begs for an exploration of new ways to perform multiloop calculations. Finally, it would be very interesting to promote our calculation to the fully non-abelian case and consider three-loop scattering amplitudes for the gg → gg process. We look forward to pursuing these lines of investigation in the future.