Two-loop master integrals for the planar QCD massive corrections to di-photon and di-jet hadro-production

We present the analytic calculation of the Master Integrals necessary to compute the planar massive QCD corrections to Di-photon (and Di-jet) production at hadron colliders. The master integrals are evaluated by means of the differential equations method and expressed in terms of multiple polylogarithms and one- or two-fold integrals of polylogarithms and irrational functions, up to transcendentality four.


Introduction
In this paper, we consider the calculation of the two-loop Master Integrals (MIs) needed for the evaluation of the massive NNLO QCD planar corrections to the hadro-production of a photon pair (di-photon production). The same two-loop masters enter the evaluation of the massive NNLO QCD planar corrections to the hadro-production of two jets (a pair of gluons, qq or q(q)g pairs, di-jet production). For massive QCD corrections we mean QCD corrections that contain a loop of heavy-quarks.
The two processes are very relevant for the physics programme at the LHC. Di-photon production, due to an experimentally clean final state, provides an important test of the Standard Model (SM) and constitutes an irreducible background for a Higgs, produced in gluon fusion, that decays into two photons. The current theoretical description of the di-photon production includes the NNLO QCD corrections [1] due to massless states, and a part of the N 3 LO massless corrections in the gluon-gluon channel [2], while the massive corrections (at the two-loop level) are not included. 1 Di-jet production is a dominant process at the LHC. It enters the determination of the strong coupling constant, it is sensible to the value of parton distribution functions and it is important for searches of new physics, beyond the SM. Very recently the NNLO QCD corrections in the purely gluonic channel were computed [4]. In [5], a phenomenological study of the dijet production, doubly differential in the mass of the dijet system and in the rapidity difference, was presented. The study involves NNLO QCD corrections to all the partonic channels, leading in the JHEP01(2018)048 number of colors. These perturbative corrections include massless states and the massive corrections (at the two-loop level) are not considered in the analysis.
If we expand the di-photon production cross section in powers of the couplings, the fine structure constant α and the strong coupling constant α S , the order α 2 α 2 S is the first perturbative order at which QCD massive corrections appear. In the qq channel, this gives a genuine two-loop correction, interference between the two-loop and the tree-level matrix elements. In the gg channel, since the photons do not couple to the gluons, the order α 2 α 2 S contains only the interference between one-loop diagrams, in which a heavy-quark box mediates the coupling between gluons and photons. Two-loop corrections appear, in this channel, at the following perturbative order, α 2 α 3 S and arise from the interference of two-loop with one-loop diagrams.
For the di-jet cross section, the leading order (LO), O(α 2 S ), is present in both the partonic channels. Therefore, NNLO QCD corrections contain genuine two-loop corrections initiated by a qq, gg or q(q)g pairs. In the case of the partonic process qq → qq the massive corrections involve a closed heavy-quark loop as correction to the gluon propagator. For the partonic processes qq → gg, gg → qq and q(q)g → q(q)g, the massive corrections involve instead actual two-loop massive box diagrams, as the ones used for the di-photon production. The channel gg → gg needs additional massive box diagrams, with respect to the ones presented in this paper; some of them are known already in the literature [6].
The di-photon production partonic cross section has a simple color structure. At the NNLO, it is proportional to C F C A , where C F is the Casimir of the fundamental representation of SU(N c ) and C A is the Casimir of the adjoint representation. This implies that both planar and crossed Feynman diagrams contribute to the sole gauge-invariant color coefficient. The di-jet production partonic cross section has, instead, a complicated structure in terms of color coefficients. In this case, the planar diagrams contribute to all the color coefficients, while the crossed diagrams do not enter the leading color one. This means that, the computation of the planar diagrams alone allows, nevertheless, to give physical predictions in the leading color approximation, while for the di-photon case, this is not possible.
We approach the calculation using Feynman diagrams. We construct the interference between the two-loop and the tree-level amplitudes; this is expressed as a combination of dimensionally regularized scalar integrals. These integrals are reduced to a set of 37 MIs using the computer programs 2 FIRE [11][12][13] and Reduze 2 [14,15], that implement integration-by-parts identities [16][17][18] and Lorentz-invariance identities [19].
The MIs are computed using the differential equations method [19][20][21][22][23]. The system of differential equations obeyed by the MIs is cast in canonical form [24] (see also [25][26][27][28][29][30][31][32][33][34]). The solution is expressed in terms of Chen's iterated integrals, represented whenever possible in terms of Goncharov's polylogarithms (GPLs) [35][36][37]. The part at transcendentality four of seven four-point functions and one part at transcendentality three are given in terms of single and double parametric integrations. The reason is that it was not possible to find a change of kinematic variables able to "linearize" the whole set of square roots that belong to the original alphabet.

JHEP01(2018)048
The analytic results presented in this article are collected in ancillary files which we upload with the arXiv submission. 3 The article is structured as follows. In section 2, we give our notations and conventions. In section 3, we discuss the system of differential equations and the canonical form of the set of MIs. Moreover, we present the alphabet that enters the solutions of the differential equations and the transformation of variables that linearizes all the square roots except one. In section 4, we present our basis and the transformation matrix between the MIs in canonical and in "pre-canonical" form. Finally, in section 5 we conclude. In appendix A, we discuss the prescription for the linearization of the square roots and in appendix B we give the routing of the MIs in pre-canonical form.

Notations
We consider the basic processes 4 qq → γγ, (gg, qq) and the crossed q(q)g → q(q)g, in which the initial partons have momenta p 1 and p 2 and the final photons or partons have momenta p 3 and p 4 . The external particles are on their mass-shell p 2 i = 0. We introduce the Mandelstam variables such that s + t+ u = 0. Since we consider 2 → 2 scattering processes with massless external particles, the physical region is defined through the following relations where θ, 0 < θ < π, is the scattering angle in the partonic center of mass frame. For later convenience we define the following dimensionless ratios where m t is the mass of the heavy quark that runs into the loops. The NNLO QCD planar corrections to the partonic processes listed above can be calculated reducing to the MIs two topologies, shown in figure 1 and defined by the integrals where a i , with i = 1, . . . , 7, are integer numbers, while a 8 and a 9 are natural numbers. The D i , i = 1, . . . , 7, are the denominators involved and D 8 , D 9 are the numerators, d is the dimension of the space-time, = 4−d 2 , γ E = 0.5772.. is the Euler-Mascheroni constant and the normalization is such that The results expressed in terms of GPLs are provided in text files in GiNaC format, while the single and double parametric integrations are provided in ".mx" format, readable by Mathematica 11 [38]. 4 The processes gg → γγ, (gg, qq), involve additional topologies and additional MIs with respect to the ones we are going to present in this paper. The routings that we used for the two topologies are the following:

JHEP01(2018)048
The momenta k 1 and k 2 are the loop momenta. The number of MIs for the two topologies is 32 for Topology A and 25 for Topology B. However, some of them appear in both topologies. Therefore, the total number of independent MIs is 37. Among them, 12 are new and presented in this paper for the first time.
Topologies A and B, do not complete the possible integrals involved in the planar QCD corrections to qq → γγ. Among them, there are three-point functions that are not included in the subtopologies of topologies A and B. These are the three-point functions that occur for instance in the calculation of the NLO QCD corrections to the Higgs production in gluon fusion (or Higgs decay into a pair of photons) and the corresponding MIs were studied in [39,40].
If we denote by f ( x, ) the vector of MIs found after the reduction process, we find that f ( x, ) satisfies a system of first order linear differential equations with respect to the kinematic invariants x, that we can write in differential form: In eq. (3.1), A( x, ) is the matrix associated with the system, which in general depends on the kinematic variables x and on the dimensional regulator . The matrix A( x, ) satisfies the integrability conditions

JHEP01(2018)048
It was argued in [24], that it is possible to choose a particular basis of master integrals, the so called Canonical Basis, in which the system of differential equations can be cast into the following simplified form: dÃ( x) is a d-log one-form, i.e. the entries ofÃ( x) are Q-linear combinations of logarithms.
In the case in which the canonical basis can be found, the master integrals can be written in terms of multiple polylogarithms.
In the canonical basis the solution of the system is formally given in terms of Chen iterated integrals [41] f where P stands for the path-ordered integration, γ is some path in the kinematic invariants space and f 0 ( ) is a vector of boundary conditions. We distinguish among two cases: the first one is the case in which the matrixÃ( x) is rational in the kinematic invariants, while in the second case the dependence on the kinematic invariants can be algebraic, i.e.Ã( x) depends also on roots of x.
In the first scenario, the entries are linear combinations of terms of the kind log(x i −α k ), where α k are algebraic functions of kinematic invariants. In this case, once specified a path γ := γ(t), the solution can be written explicitly in terms of Goncharov's multiple polylogarithms [35,36] In the second scenario, the entries contain roots of algebraic functions of kinematic invariants. In the case in which all the roots can be linearized simultaneously with a change of variables, the solution can be found, again, in terms of Goncharov's multiple polylogarithms. If it is not possible to find such a transformation, one is left with repeated integrations over rational functions and/or multiple polylogarithms and a direct integration in terms of multiple polylogarithms is not easy. In this case, a possible strategy to arrive at the solution is to use the concept of symbol of an iterated integral in order to construct an ansatz in terms of multiple polylogarithms of a certain weight (see for instance section 3.2 of [42]). Another strategy, is to follow the idea outlined in [6]. If the weight-2 is known in terms of Goncharov's multiple polylogarithms, we can integrate the following two integrations (up to weight 4) numerically. Moreover, integrating by-parts, we can reduce the double integration to a single integration, that involves irrational functions times weight-3 polylogarithms. For the process under study, we found a transformation that linearizes a part of the roots present in the system of differential equations, but not all of them. In the set of new variables, the weight-2 is expressed in terms of multiple polilogarithms. It is therefore possible to follow the approach just outlined and to find the solution in terms of one-and two-fold numeric integrations.

The alphabet
The matricesÃ( x) of the systems of differential equations for the two topologies A and B depend on the following two groups of five square roots: Regarding Topology B, we managed to further split system B into two subsystems B 1 and B 2 that do not share master integrals, and have a smaller set of square roots: which we could linearize simultaneously using the transformations (see appendix A) . For what concerns Topology A, it was not possible to split system A into independent subsystems with a smaller number of roots. We managed to linearize simultaneously the first four square roots, while it was not possible to linearize the fifth. The presence of the fifth square root in eq. (3.7) can be found just in the differential equations that involve the MIs of a 5-denominator four-point topology and the four-point topology at 7 denominators. In the rest of the system, only the first four square roots appear. This makes in such a way that we can separate the matrix of the system,Ã(u, v) into two pieces (3.13) The first termÃ nr (u, v) contains the same set of roots of the subsystem B 1 , whileÃ r (u, v) contains also the additional root u(u + 8uv + 16(1 + u)v 2 ). We performed the same change of variables used for the subsystem B 1 . With this choice we could linearize the set of square roots ofÃ nr (u, v) while the termÃ r (u, v) still contains square roots in the new kinematic variables w and z. This splitting of the system is functional to the solution. Indeed, the MIs whose differential equations are described by the entries of the matrix A nr (u, v) are analytically expressible in terms of GPLs, to all orders in the expansion. On

JHEP01(2018)048
the other hand, the MIs whose differential equations are described by the matrixÃ r (u, v) are not analytically expressible in terms of GPLs and involve one-and two-fold numeric integrations.

The master integrals
In order to find the canonical basis for the MIs several approaches exist [24,25,29,[32][33][34]43]. In this paper we adopt the semi-algorithmic method described in [28], which is based directly on the analysis of the system of differential equations. We put the system into canonical form using a "bottom-up" approach; we start from the subtologies with lowest non-trivial number of denominators and we proceed with subtopologies

JHEP01(2018)048
with higher number of denominators, up to the completion of the rotation of the whole system from pre-canonical to canonical form. In general, we proceed following three steps: • Step 1. If we denote with f i (x) the MI under consideration, we choose the powers a i of the denominators in order to get the differential equation that concerns f i (x) in the following form The dependence of the non-homogeneous term Ω ij (x, ) has to be of the form where p a are real numbers.
• Step 2. We remove the term H 0,ij (x) rescaling the master integral f i (x) by a functions h 0,ij (x) which satisfies the differential equation where, again, the indices run over coupled master integrals. For all the master integrals that we studied we found that the functions h 0,ij (x) are algebraic in the kinematic invariants. Therefore, the differential equation in the new master integral Step 3. We put in canonical form the non-homogeneous part shifting the MIf i (x) as In order to remove the first and third term in (4.2), the functionsω 0,ij (x) andω a,ij (x) must satisfy the following system of differential equations: where G 1,ij (x) is the matrix of the system of differential equations for the subtopology g j (x), which is already in caonical form: Given the set of 37 MIs in pre-canonical form, as shown in figure 2, we define the canonical basis, f 1 , . . . , f 37 , from the following relations: 14) (4.52) The initial conditions f 0 ( ) for the MIs are fixed in the point s = t = 0, where all the canonical master integrals vanish except for the double tadpole, which is known analytically, and master integrals f 2 and f 3 which are also known analytically, since they are product of two one-loop integrals.
The master integrals f 1 -f 21 and f 24 -f 26 were already present in the literature [44][45][46][47][48][49]. The master integrals f 29 , f 30 , f 31 were considered in [6]. The GPLs are found integrating the master integrals along the integration contour given by γ = γ w + γ z , where The formal integration in w and z is made in the Euclidean region 5 The formulas obtained, however, are valid in both the Euclidean and Minkowski regions, after analytical continuation, which can be done adding a small imaginary part to the Mandelstam variable s according to the Feynman prescription s + i0 + . In order to study the numerics of our analytic results in the Euclidean and Minkowski regions, we have to invert the relations (3.11), (3.12) finding w as a function of u and z as a function of v, for fixed values of u. w(u) and z(u, v) have four different branches. All the branches of w(u) are discontinuous in the points u = 4 5 , u = 1 .

JHEP01(2018)048
Therefore, in order to have the correct value of w for u spanning the complete domain, we must choose different branches of w(u) and link them continuously in the singular points (4.56). When u ∈ 0, 4 5 , w is given by and it varies in the real range w ∈ [0, 1]. When u ∈ 4 5 , ∞ , w is given by The domain u ∈ 4 5 , 1 is mapped into the real range w ∈ 1, 1 + √ 2 , while for u > 1 the variable w becomes complex. Regarding the variable z, we find that for any fixed value of u we can map the whole Euclidean domain for v into a closed finite range for z using a single branch of z(u, v). The analytical continuation to the Minkowski region (u < 0) is performed, according to the Feynman prescription, adding a small imaginary part to the Mandelstam variable s, s + i0 + . In particular, the range u ∈ [0, −1] is mapped into a complex region for w, using the expression (4.57). Similarly, in the range u ∈ [−1, −∞), using the expression (4.60), we find that u is mapped into a complex region of w.
We performed numerical checks of our GPL analytic results, comparing them with the results obtained using the software FIESTA4 [50][51][52] for points both in the Euclidean and Minkowski regions, finding complete agreement.
As mentioned previously, the appearance in the alphabet of the Topology A of the fifth root u(u + 8uv + 16(1 + u)v 2 ) makes not possible to linearize simultaneously all the roots of the system. Therefore, we splitted the differential equation matrix into a sum of a part containing only GPL terms and a part that contains the root u(u + 8uv + 16(1 + u)v 2 ). The numerical integration is performed exploiting the fact that the weight 2 can be written in terms of GPLs for the whole system. This is because the fifth root affects the solution only from weight 3 on. Hence, we wrote the weights 3, f i (x) and 4, f (4) i , as single and double numerical integrations, respectively, over an analytic kernel of weight 2 For the numerical integrations, we chose the same paths that brought to the analytic expressions in terms of GPLs. We note that, in so doing, the integration in w gives rise to JHEP01(2018)048 expressions that are cast in polylogarithmic form. The sole numeric integration remains the one in z. We checked our numerical results, coming from the numeric integrations, against the software FIESTA4 for points in the Euclidean region, finding complete agreement. The analytical continuation to the Minkowski region of the one-and two-fold parametric integrations is not considered in this article.

Conclusions
In this article we presented the calculation of the master integrals needed for the evaluation of the NNLO QCD planar corrections to di-photon (and di-jet) production in hadronic collisions.
The system of differential equations satisfied by the master integrals is cast in canonical form and solved in terms of Chen's iterated integrals. We represented these integrals in terms of Goncharov's multiple polylogarithms (GPLs), whenever it was possible to linearize the set of square roots appearing in the corresponding alphabet. However, the finite parts of seven four-point functions, and a single pole of one of them, still involve one-and twofold numeric integrations of polylogarithms multiplied by irrational functions. Choosing as a path for the iterated integration the same path that brought to the representation in terms of GPLs, we were able to integrate the letters that are function of the variable related to the partonic c.m. energy in closed form in terms of GPLs, remaining with an integration in the sole variable related to the partonic momentum transfert.
The part expressed in terms of GPLs is integrated formally in a small part of the Euclidean region. However, it can be evaluated in the whole Euclidean region and analitically continued to the Minkowskian region simply adding the correct causal prescription to the Mandelstam invariant s. We checked that the numerical values obtained from our analytic expressions using the routines in [53] agree with the ones obtaines with FIESTA4.
Concerning the numeric integrations, the fact that they involve only the variable related to the momentum transfert, allows for a strightforward evaluation in the whole Euclidean region and for a possibly simple analytical continuation to the Minkowski region, that, nevertheless, we do not consider in the present work.
With the master integrals presented in this article (and other already present in the literature), it is possible to evaluate the leading color contribution of the massive corrections to the di-jet production. It is not possible to evaluate a gauge independent quantity in the di-photon production. For this we need the contribution of the crossed diagrams, that will be considered in a subsequent publication.

A Linearization of square roots
In order to write a solution in terms of multiple polylogarithms, the dependence ofÃ( x) on the kinematical invariants has to be rational. Therefore, we need to find a change of variables that linearizes simultaneously all the roots in the matricesÃ( x). In this appendix, we describe the method that we used in order to obtain the transformations (3.11), (3.12), needed for the complete linearization of square roots belonging to Topology B.
Once the transformations that linearize the set of roots (A.1) are found, their inverse linearizes the set of roots (3.8).
The problem of roots linearization is connected to the diophantine equation, which is a polynomial equation in many variables such that only its integers solutions are studied. It is possible to find a change of variables that linearizes simultaneously two or more roots solving iteratively the diophantine equation associated to them. The idea is to transform the argument of the root into a perfect square, exploiting the solution of the diophantine equation associated as a parametrization for such transformation.
Consider the roots The second solution is exactly the transformation ofũ that linearizes the root √ 1 +ũ, as can be easily checked.
Once the first root is linearized, we consider the second root, √ũ − 1. We transform the original variables according to the transformation that linearizes the first root, eq. (A.6). The diophantine equation associated to the root √ũ − 1 is