Two-loop form factors for diphoton production in quark annihilation channel with heavy quark mass dependence

We present the computation of the two-loop form factors for diphoton production in the quark annihilation channel. These quantities are relevant for the NNLO QCD corrections to diphoton production at LHC recently presented in arXiv:2308.10885. The computation is performed retaining full dependence on the mass of the heavy quark in the loops. The master integrals are evaluated by means of differential equations which are solved exploiting the generalised power series technique.

More recently, scattering amplitudes belonging to higher orders in the strong coupling α s (i.e.beyond the NNLO) have become available (in the massless case): the three-loop matrix element [37]; the two-loop scattering amplitudes for a photon pair in association with one jet in the leading colour approximation [38][39][40], and, very recently, the full colour case [41].The two-loop scattering amplitudes for diphoton production in gluon fusion [42] together with the recent computation of diphoton production in association with one jet [43] at NNLO emphasise that all the building blocks are in place for the next-to-next-to-next-to-leading order (N 3 LO) massless calculation.However, the implementation of slicing subtraction methods to reach the N 3 LO accuracy could be challenging in this case, due to the presence of a high number of particles in the final state and the use of a photon isolation prescription.More clearly, the NNLO calculation of the diphoton cross section in association with one jet at small diphoton transverse momentum could be very CPU demanding.
In the massive case, the first non-trivial corrections appear at the NNLO, through the inclusion of top quark loops and top quark radiation 1 .The mass effects of the so-called box contribution (gg → γγ) were discussed in ref. [23] (together with partial N 3 LO contributions).Due to the large gluon luminosity at the LHC, the size of the box contribution is of the order of the Born subprocess q q → γγ.It is therefore of interest to calculate the corrections of the following perturbative order (i.e.N 3 LO contributions).Regarding the gluon fusion channel only, the simplest approach (which captures very sizable components) is to consider the NLO QCD corrections to the box contribution, since these form a gauge invariant subset [42] of the whole N 3 LO gluon fusion channel.In this context, two recent papers have shown the impact of these massive NLO QCD corrections on the gluon fusion channel [44,45].
Considering the full NNLO accuracy, there are still three missing ingredients that were not available or not presented together in previous phenomenological studies in the literature: i) the massive one-loop real-virtual contribution (q q → γγg and q(q)g → γγq(q)), the double real radiation of top quarks (q q → γγt t and gg → γγt t) and the two-loop virtual corrections to the Born subprocess q q → γγ.
In this paper we report on the computation of the two-loop form factors for diphoton production in the quark annihilation channel, where the full dependence on the heavy quark mass, which appear in the loops, is retained.This contribution is included in the recently presented phenomenological study of full massive NNLO QCD corrections to diphoton production [1].
The two-loop form factors have been computed employing standard techniques for scattering amplitude calculations.We considered the partonic process q q → γγ and the relative two-loop Feynman diagrams, which contain a massive heavy quark loop, as shown in figure 1.The associated amplitude is decomposed into a combination of tensors multiplied by scalar form factors, as described in [46], and the scalar integrals appearing in the expressions of the form factors are written in terms of a basis of master integrals (MIs).The decomposition in terms of MIs is performed using Identity-by-Parts (IBPs) relations [47,48], via the Laporta Algorithm [49], implemented in the computer code2 KIRA [57,58].
The MIs relevant for this process have been computed by means of the differential equations method [59][60][61][62][63][64][65][66].While the integrals associated to the planar topologies have already been studied in the literature [67][68][69][70][71][72][73], a complete computation for the non-planar ones is still not available.Indeed, while the planar MIs admit an analytic solution in terms of Multiple Polylogarithmic functions (MPLs) [74][75][76][77][78][79], it is known that the functional space for the analytic solution of the non-planar double-box family [80] contains elliptic integrals [81][82][83][84][85][86][87][88][89][90][91][92].In recent years, a big effort has been devoted to the understanding of the analytic structure of Feynman integrals which do not admit an expression in terms of MPLs.However, even in the cases in which an analytic solution in closed form is available, the numerical evaluation of the functions associated to such solution can be extremely challenging for phenomenological applications [93,94].
In order to be able to overcome the issues previously described, we choose to exploit the generalised power series method [95][96][97][98][99][100][101] to solve the systems of differential equations associated to the MIs.This technique, currently implemented in two public computer codes [102,103], has recently attracted a lot of interest due to its wide range of applicability and it has been successfully employed in several phenomenological applications [1,[104][105][106].Specifically, in the calculation reported in this paper, we used the software DiffExp [102] to obtain a semi-analytic solution for the MIs.
The paper is organised as follows.In section 2, we describe the general setup for the computation and we set the context in which the two-loop form factors, presented in this paper, are relevant.We discuss the form factors computation along with their UV singularity structure and the renormalisation procedure.In section 3 instead we report on the MIs calculation.We describe the different integral families that appear in the computation and the approach that we used to solve the differential equations, along with a brief analysis on the geometry underlying the actual analytic solution.
Moreover, as ancillary material attached to this paper, we furnish the analytic expressions of the finite reminder for the form factors, alongside with Mathematica files which allows for a standalone evaluation of the MIs with DiffExp.

Computational setup and amplitude structure
In this paper, we consider the two-loop form factors for diphoton production in the quark annihilation channel with a heavy quark loop.At the partonic level, the scattering amplitude proceeds as the Born subprocess: q(p 1 ) + q(p 2 ) → γ(p 3 ) + γ(p 4 ). (2.1) The kinematics for this process is described by the Mandelstam variables3 where the external particles are on-shell, i.e. p 2 i = 0, and we indicate with m 2 t the heavy-quark squared mass 4 .In order to obtain the scattering amplitude, we generated the relevant Feynman diagrams using the FeynArts package [108].We found a total number of 14 diagrams contributing to the amplitude, the representative ones are shown in fig. 1.We write the scattering amplitude in terms of form factors, which are decomposed into a basis of 72 MIs exploiting IBPs reduction [47-49, 51-57, 109], as implemented in the software Kira [57].
The MIs contributing to this process can be described by three different scalar integral topologies (modulo exchange of the two final photons).Specifically, the MIs for the Feynman diagrams (a) and (b) in fig. 1 are associated to the integral families PLA and NPL, respectively, as defined in section 3. Similarly, the MIs for the diagrams (c), (d) and (e) can be grouped into one scalar integral family, PLB, also defined in section 3. The MIs of the families PLA and PLB were already known in the literature [67].Regarding the non-planar topology NPL, while most of the MIs have already been studied [67,72,73,97,[110][111][112], the double-box top-sector have not been considered in the literature yet, and therefore its computation represents an original result by itself.
For this project, we performed an independent calculation of all the MIs by means of the differential equations method [59][60][61][62][63][64][65][66].In particular we solved the system of differential equations semi-analytically exploiting the generalised power series expansion technique, as described in [100] and implemented in the software DiffExp [102].
The two-loop amplitudes computed in this paper constitute a necessary ingredient of the recently presented full massive NNLO QCD corrections to diphoton production at hadron colliders [1].We anticipate here, that our two-loop form factors, taking into account the full dependence on the top quark mass, are finite (free of IR divergences after UV renormalization).Therefore, they can be included directly (without the use of any IR regularization prescription) in any numerical implementation of the NNLO cross section.In the following paragraph we will illustrate the precedent statement with a specific example based on the q T -subtraction method [113,114] (which can easily be extended to any other subtraction method).
To this end, we consider the following scattering process, where h 1 and h 2 are the colliding hadrons.At the NNLO, the cross-section for this process can be computed using the q T -subtraction method [113][114][115] as follows (2.4) The terms inside the square brackets dσ γγ+jets NLO and dσ CT NLO represent the cross section for diphoton plus jet production at NLO [36] and the corresponding counterterm, needed to cancel the associated singularities in the small-q T limit.The coefficient function H γγ NNLO (defined in [115]) is the so-called hard-virtual function and it includes the one-loop and two-loop corrections to the Born subprocess.This object admits a perturbative expansion in terms of the strong coupling α S : In our particular case (diphoton production) the one-loop contribution to eq. (2.5) was calculated in [33], while the massless two-loop contribution was first calculated in [35] and later in [37].The explicit expressions of the hard virtual factor (computed with the previous massless oneand two-loop amplitudes) in the hard resummation scheme are given in Appendix A of ref. [114].The inclusion of the new two-loop massive form factors proceeds by simple addition to that of the massless case [35].After regularising the IR divergences (e.g. in the in the hard resummation scheme) [114] present in the massless two-loop amplitude [35], our massive two-loop contribution can be added directly to the massless finite remainder.
After generating all the relevant Feynman diagrams for the process, the amplitude A q q,γγ has been decomposed in terms of form factors [46], and the UV singularities have been regularised in dimensional regularisation.The expression obtained has been used to compute the NNLO corrections to the hard function H γγ coming from two-loop diagrams which involve a massive top quark loop.Specifically, from the knowledge of the finite remainder A (fin) q q,γγ of the amplitude, we can obtain the hard function from the all-orders relation [114]: where q q,γγ is the Born-level amplitude for this process: with α em , the QED coupling, Q q is the electric charge of the incoming quarks and N c is the number of colours.We also performed a sum over initial and final polarisations and initial colours.This massive hard-virtual coefficient represents the last missing ingredient necessary to perform a NNLO phenomenological study, for diphoton production at the LHC, which takes into account the complete dependence on the top quark mass [1].

Form factors
The bare scattering amplitude can be written 5 as where δ ij is the Kronecker delta function with i and j the color indices of the incoming q q pair, ϵ * µ λ3 (p 3 ) and ϵ * ν λ4 (p 4 ) are external photon polarisation vectors and v s2 (p 2 ), u s1 (p 1 ) the quark spinors.The amplitude (2.8) can be further decomposed in terms of a set of four independent tensors 6[46] which are built using external momenta and polarisation vectors: where the T k are chosen as: (2.10) The decomposition (2.9) has been achieved also by enforcing the physical conditions ϵ λi • p i = 0 for the polarisation vectors, and by choosing as reference vectors for the external photons the condition ϵ λ3 • p 2 = ϵ λ4 • p 1 = 0, which implies for the photons the polarization sums: where g µν = diag(1, 1, 1, 1).The coefficients of T k are the so-called scalar form factors, F k .These objects are functions of the kinematic invariants of the process, and of the space-time dimension D, and they can be written in terms of scalar Feynman integrals.Their expression can be obtained by applying a set of projectors, {P k }, k = 1, 2, 3, 4, to the amplitude A q q,γγ : where (2.13) The form factors F k admit the following perturbative expansion: where α S is the strong coupling constant.At leading order we have: (2.15) The massive corrections, we are interested in, appear starting from the two-loop order.Therefore, they affect only the term F k .For the rest of this paper we will focus just on this contribution, which we will refer to as F (2) k,top .We generated all the relevant Feynman diagrams for this contribution using FeynArts [108], and we found 14 different two-loop Feynman diagrams, which can be grouped into five different categories, as depicted in fig. 1.We used FORM [116,117] to apply the projectors to the Feynman diagrams and perform the Dirac algebra.The form factors, then, were expressed as a linear combinations of 72 MIs, {J i }, which are defined in Appendix A.
We find the following expression for the form factors F (2) k,top : where 2Nc is the Casimir of the fundamental representation of SU (N c ) and Q t is the electric charge of the top quark running in the loop.The two contributions F We perform our computations in the context of dimensional regularisation.As a consequence, potential ultraviolet (UV) and infrared (IR) singularities can appear in the form factors as poles in the dimensional regulator ϵ = (4 − D)/2.However, since the diagrams with a top loop start contributing to the q q channel at the two-loop order, F k,top does not have IR singularities and therefore all the ϵ poles are of UV origin.Furthermore, F k,top;2 is also free of UV divergences and then the UV poles come only from the contribution F (2) k,top;0 .

UV Renormalisation
We renormalise the bare form factors in a mixed scheme.The external quark fields are renormalised on shell; the strong coupling constant is renormalised in a scheme in which the light-quark contribution is treated in MS, while the heavy-quark contribution is renormalised at zero momentum.
No renormalisation is needed for the top quark mass m 2 t at this order in perturbation theory and the same occurs for the external photon field.
We have, then where the superscript "R" stands for "renormalised".The renormalisation factors Z q and Z α S admit a perturbative expansion in the strong coupling constant [118][119][120][121]: (2.18) In our renormalisation scheme, we have α,N l ,MS + δZ where µ is the renormalisation scale, N l is the number of light flavours, in this case N l = 5, N h is the number of heavy quarks, in this case N h = 1, the quantity C A = N c is the Casimir of the adjoint representation of SU (N c ), and T F = 1 2 .Therefore, the renormalised form factors F (2) R k,top are determined by the following formula: k,top + δZ (2)  q F (0) k . (2.21)

Master Integrals Computation
In this section we discuss the details of the MIs computation.Specifically, we define the three scalar integral topologies which describe all the MIs that appear in the amplitude (modulo exchange of final photons).The MIs are computed through the differential equations method.The system of differential equations associated to the MIs is solved semi-analytically employing the generalised power series expansion technique, as described in [100] and implemented in the software DiffExp [102].Finally, we show that some of the MIs that appear in the computation admit an analytic solution in terms of elliptic functions.
The system of differential equations for the planar families, PLA and PLB, is in canonical form, while for the non-planar family, NPL, we have canonical differential equations just for the sectors whose analytic solution could be given in terms of MPLs.We notice that the bases of MIs in which we solve the differential equations, which we will refer in the following as ⃗ f for the planar families and ⃗ g for the non-planar one, and the one in which we write the form factors are different.We choose this approach in order to avoid square roots of the kinematic invariants in the expressions of the form factors.In the ancillary material attached to this paper, we furnish the rotation matrices from the bases ⃗ f and ⃗ g to the basis {J i } defined in Appendix A.

Denominator
Integral family PLA Integral family PLB Integral family NPL Table 1: Routing definition for the three scalar integrals families PLA, PLB and NPL.

General Setup
The three integral families which describe the relevant MIs for this process are defined as follows: where topo ∈ {PLA, PLB, NPL} labels the families and the scalar products D 1 , ..., D 9 are given in table 1.The indices n 1 , • • • , n 7 are non-negative while the indices n 8 and n 9 are non-positive.The computation is done in dimensional regularization with D = 4 − 2ϵ dimensions, and our convention for the integration measure is where µ is the dimensional regularization scale.
With the definitions given in table 1, all the scalar Feynman integrals appearing in the amplitude can be mapped to one of the three integral families, or families obtained by permutation of the external momenta.The MIs appearing in the form factors are listed in appendix A and shown in fig. 2.

Planar Families PLA and PLB
The scalar integrals family PLA is described by a set of 32 master integrals ⃗ f (⃗ x, ϵ) 7 , where is the vector of the kinematic invariants with respect to which we derived the differential equations.
The system of differential equations for these integrals has been derived in canonical logarithmic form [65] in ref. [67]: where d is the total differential with respect to the kinematic invariants, and the matrix A(⃗ x) is written as linear combinations of logarithms: The c i represent matrices of rational numbers, while W P LA = {w i (⃗ x)} is the alphabet of the solution and it is made by algebraic functions of the kinematic invariants.Specifically, the alphabet is made by the following set of 21 letters: where r 1 , • • • , r 5 are a set of square roots of the kinematic invariants: (3.7) The knowledge of the logarithmic canonical form of the differential equations (3.4), together with the alphabet structure (3.6) would allow us, in principle, to obtain a fully analytic representation for the system of the MIs.However, the presence of the set of square roots given in Eq. (3.7) makes the achievement of such analytic expression non trivial.Indeed, these square roots are not simultaneously rationalizable.As a consequence, in order to obtain a fully analytic representation of the solution one would have to exploit symbol level techniques [122,123].For the purpose of this project, we found that the semi-analytic evaluation, which we achieved exploiting a generalised power series expansion method, was sufficient to perform phenomenological studies.
The boundary conditions for the system are provided in the origin of kinematic variables y = z = 0, where all the MIs vanish except for the two masters f 1 and f 2 , for which we use the following analytical expressions: where J 1 and J 2 are the pre-canonical masters shown in fig. 2 and defined in Appendix A.
Regarding the second planar family, PLB, we observe that we do not need to set up a system of differential equations for it.Indeed, every master integral of this family, except J 21 , is equal to one of the MIs of the two other families (modulus permutations).For J 21 , by integrating analytically its differential equation, we obtain the following analytical expression: This formula was cross checked analytically with [124] and numerically with AMFlow [125].

Non-Planar Family NPL
The non-planar scalar integrals family, NPL, is more complicated to study with respect to the two planar families.The total number of MIs ⃗ g8 associated to this topology is 36.The system of differential equations can be divided in two different subsets: • (I) Canonical logarithmic: the subset of MIs whose differential equations is in canonical logarithmic form; • (II) Elliptic sectors: this subset contains MIs whose analytic solution involve elliptic functions and it is not written in ϵ-factorized form.
The system of differential equations for the subset (I) has been put in canonical logarithmic form.
We choose for them the same normalization as in [112]: The second sector whose analytic structure is characterised by the presence of elliptic functions is the top sectors of the topology NPL, i.e. the double-box integral, shown in fig. 3 (b), which contains 4 MIs: g 33 , g 34 , g 35 and g 36 .We choose, as basis for this sector, the following scalar integrals: With this choice of normalization, initial conditions in the origin y = z = 0 are: where f 1 and f 2 are given in Eq. (3.8).
The non-polylogarithmic structure of this sector is two-fold.First, the differential equations for the MIs g 33 , g 34 , g 35 and g 36 contain the triangle integrals g 31 and g 32 in the non-homogeneous part of the system.As a consequence the analytic solution of the differential equations requires the integration over kernels that contains eMPLs.Moreover, also the homogeneous part of the differential equations itself contains elliptic functions.This statement can be verified by studying the maximal cut of the double-box integral [66]: (3.15) In order to perform such computation it is convenient to follow the loop-by-loop analysis [127,128] of the Baikov [129] representation of the integral I db .For reader's convenience, we briefly introduce the basic idea of Baikov representation and its application to maximal cut analysis, and we refer to refs.[127][128][129][130] for a proper description of the method.The Baikov representation of an L-loop Feynman integral in D dimensions, with E independent external momenta, is written as: The variables z i are just the scalar products D i which define the integral topology.Indeed, the basic idea of the Baikov representation is to write a scalar Feynman integral in a form where the integration variables are exactly the scalar products associated to the topology.The factors U and P are, respectively, the Gram determinant of the external momenta p 1 , • • • , p E : and the Gram determinant of the external momenta together with the loop ones: where G ij (⃗ v) = v i • v j is the Gram matrix.The factor C L E is a normalization constant, which is immaterial for the present discussion.
In the Baikov representation the maximal cut of the integral I(a 1 , • • • , a n ) can be calculated from the multivariate residues of the expression (3.16) [131].This operation is more efficiently done within the so-called loop-by-loop approach [127].In a nutshell, in this approach we split the integral under study in subsequent one-loop integrals and then we take the maximal cut of each sub-integral sequentially.Let us consider explicitly our double box integral (3.15).By inspecting the integrand of (3.15) we see that it can be split into two one-loop box integrals as follows: where z ± = 1 2 −s − 2t ± √ s s + 16m 2 t .As we can see from Eq. (3.20) the result of the maximal cut for the double box non planar integral in fig. 3 (b) is a one-fold integral.It is possible to show that Eq. (3.20) can be written in terms of complete elliptic integrals of first and third kind, where the elliptic curve is given by the polynomial of fourth-order in the integration variable z 8 : Remarkably, we observe that the elliptic curve (3.21) degenerates to the same curve of the massive triangle in fig. 3 (a) [93,112,126] in the forward limit t = 0.

Semi-Analytic solution with Generalised power series
We conclude this section by describing the solution for the systems of differential equations associated to the planar topology PLA and the non-planar topology NPL.As already mentioned, we choose to exploit the generalised power series method, as described in [100] and implemented in the software DiffExp [102], to obtain a semi-analytic solution for the set of MIs.This method has the advantage of not being limited by the functional space in which the MIs would be analytically represented.This feature allows us to avoid the issues connected with the presence of MIs which admit an analytic solution in terms of elliptic integrals, for which both the understanding of their analytic structure and the numerical evaluation can still represent a bottleneck for phenomenological applications.
We exploit the method to build a grid of points for the contribution of the corrections considered in this paper to the hard function H γγ N N LO .After interpolation, the grid has been used in the fully massive NNLO phenomenological study for diphoton production in [1].The grid has been generated directly in the physical region of the phase-space for this process: where θ, 0 < θ < π, is the scattering angle in the partonic center of mass frame.Since the evaluation time, within DiffExp, for the MIs needed in this process is relatively low, we can build the grid of points as follows 9 .We consider a total number of 13752 points in the following range for the scattering angle θ and the energy of the center of mass s: The points p i,j of the grid are equally spaced in the coordinate system (s, t) 10 as follows: (3.24) 9 We notice that for more CPU demanding computations more refined approaches have to be used. 10For the purpose of this discussion we use the dimensional Mandelstam variables defined in Eq. 2.2 From this point, at fixed value of s, we move along the θ axis, from cos θ = −0.99 to cos θ = 0.99, up to the point p 0,23 .Then, we increase the value of s and we move in the other direction along the θ axis up to the point p 1,0 , and so on so forth.In order to optimise the grid generation, for each evaluation we use as boundary conditions the value of the MIs obtained at the previous point.This procedure effectively increases the efficiency of the evaluation for the MIs.In particular, we managed to evaluate the MIs in all the points of the grid, with a 16 digits accuracy, in 2.5 hours for the system PLA and 10.5 hours for the system NPL, on a single core laptop.Finally, in order to validate our results, we performed numerical checks for the MIs against independent numerical evaluations done with the AMFlow package [125], which implements the auxiliary mass flow method [132,133].The MIs have been checked for several points in the physical phase-space region, finding an agreement between the two independent evaluations up to 200 digits of accuracy.As a proof of concept of the numerical checks we show in figure 5 our results for the double-box MIs in the non-planar topology NPL.

Conclusions
In this paper we presented the computation of the two-loop form factors for diphoton production in the q q channel, where the full dependence on the top quark mass has been retained.This computation represents the only missing ingredient, at two loops, in order to be able to perform a phenomenological study for diphoton production at NNLO [1] which fully takes into account the dependence on the top quark mass in all the relevant channels.
The non-planar topology which contributes to this process contains two sectors of MIs whose analytic representation cannot be given in terms of MPLs.In order to be able to exploit our results for phenomenological applications, we computed the MIs by means of differential equations, exploiting the generalised power series technique.This method proves to be of great use for phenomenological applications, especially in cases where the functional space for the MIs contains not only polylogarithmic functions.

Figure 1 :
Figure 1: Representative set of two-loop diagrams with internal heavy-quark loops, which contribute at NNLO QCD corrections to diphoton production in the quark annihilation channel.Thin black lines represents light quarks, thick black lines heavy quarks, curly lines gluons and curby lines photons.
;2 are indeed related to different powers of the top electric charge Q t .The first contribution F (2) k,top;0 is associated to the diagrams (c), (d) and (e) in fig. 1 in which the top quark does not couple to the external photons, while the second contribution F (2) k,top;2 comes from the diagrams (a) and (b) where the top quark actually couples with the photon.

Figure 2 :
Figure 2: Basis of Master Integrals for the topologies PLA, PLB and NPL.In the figure, the MIs defined with a permutation of the external momenta are not displayed.

Figure 3 :
Figure 3: Diagrams representing the two elliptic sectors in the non planar topology NPL.Thin lines represent massless particles, while thick lines massive particles.

Figure 4 :
Figure 4: Schematic representation of the procedure exploited to optimise the grid construction within DiffExp.Red dots represents the points in which the MIs are evaluated and the blue dashed lines connect the sequential evaluations.

Figure 5 :
Figure 5: Numerical checks for the non-planar double-box MIs J 39 , J 40 , J 41 and J 42 .The plots show the order O(ϵ 0 ) of the masters.Blue and dashed red lines represent, respectively, the numerical values of real and imaginary part of the MIs obtained with DiffExp.Crossed dots represent numerical values obtained with AMFlow.The evaluations are performed for different values of √ s at fixed angle θ (in this case cos θ = 0.7).
19))where C is some overall normalisation constant.U 2 and P 2 are Gram determinants associated to the change of variables for a one-loop box integral with loop momentum k 2 , similarly U 1 and P 1 are related to change of variables for a box integral with loop momentum k 1 .Then, the maximal cut of the expression (3.19) is obtained by taking the residue around the simple pole z 1 = • • • z 7 = 0: MCut I db = C z1=z2=z3=0 dz 1 dz 2 dz 3 dz 8 z 1 z 2 z 3 z 8