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 masters 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.


Contents 1 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 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.
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 Fig. 1 and defined by the integrals where a i are positive integers, the D i , i = 1, ..., 9, are the denominators involved, 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 Figure 1: Seven-denominator topologies. Thin lines represent external massless particles and propagators, while thick lines represent massive propagators.
The routings that we used for the two topologies are the following: 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].

The System of Differential Equations
For the analytic computation of the MIs, we use the differential equations method [19][20][21][22][23]. 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: 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 It was pointed out in [24], that in the case in which the master integrals can be expressed in terms of multiple polylogarithms, the system 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 this 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 sect. 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) . (3.12) For the subsystem B 1 the variable u transforms according to (3.11) and the variable v to the transformation (3.12). Viceversa for the subsystem B 2 .
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 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 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 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 masters 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 masterf 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: Although in principle the differential equations (4.6,4.7) can be as difficult as the starting system for the MIs, in all the cases under study we managed to solve them and the solutions are algebraic functions of the kinematical invariants.  Figure 2: Master Integrals in pre-canonical form. Internal plain thin lines represent massless propagators, while thick lines represent the heavy-quark propagator. External plain thin lines represent massless particles on their mass-shell. Some of the masters, that depend on a single external variable, can appear both as function of s or t. Because of this fact, we marked explicitly the functional dependence of the master with "s", "t" or "s, t".
Given the set of 37 MIs in pre-canonical form, as shown in Fig. 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 masters f 2 and f 3 which are also known analytically, since they are product of two oneloop integrals. The masters f 1 -f 21 and f 24 -f 26 were already present in the literature [44][45][46][47][48][49]. The masters f 29 , f 30 , f 31 were considered in [6]. The GPLs are found integrating the masters along the integration contour given by 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 . 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 with 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 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 masters 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 two-fold 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 masters 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 diphoton production. For this we need the contribution of the crossed diagrams, that will be considered in a subsequent publication.

Acknowledgments
We would like to thank Francesco Moriello and Simone Lavacca for useful discussions about the canonical form of the master integrals and the linearization of the square roots appearing in the alphabet. Feynman diagrams are drawn with Axodraw [54]. A part of the necessary algebraic manipulations was preformed using FORM [55].

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.
Instead of applying the procedure to the original set of roots (3.8), we consider the following set √ 1 +ũ, where the kinematic variablesũ andṽ are the inverse of the variables u and v, i.e.
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.
The variable p is promoted to a function p(t 1 ,ũ) of a parameter t 1 , such that p(t 1 , 0) = −1.
The easiest choice is the following linear parametrization: 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 and, in the variable t 1 , it becomes We note that the denominator of Eq. (A.8) is already a perfect square. Therefore, for our purpose it is sufficient to study the equation We proceed as for the previous root, finding the integer solutions of (A.9), which are t 1 = 0 and p = −1 and parametrizing p as p = t 2 t 1 − 1 .  If we combine the two transformations we find the transformationũ → t 2 2 + 1 (t 2 (t 2 + 4) + 5) 4(t 2 + 1) 2 , (A.14) that linearizes simultaneously the two roots √ 1 +ũ and √ũ − 1. For more than two square roots, we can repeat the procedure until all the roots are linearized. Indeed, the restriction of the method is based on the existence of integer solutions to the diophantine equation.
The root √ 1 +ṽ can be linearized with a transformation analogous to the second solution of Eq. (A.6), with a different parameter t, t v .
(A.17) Equation (A.17) is polynomial of degree 4 in the variable t u and of degree 2 in the variable t v . Therefore, the solution of the diophantine equation is more complicated with respect to the single variable case described in Eq. (3.19). However, in this case we can find the following solution: t v = 0, t u = 0, p = −2 .
(A. 18) Eq. (A.18) leads to the transformations (3.11) and (3.12), once we consider the variables u and v, that allow for a complete linearization of all the roots in the two subsystems of Topology B. Let us consider now the roots appearing in the alphabet of Topology A. Beside the roots observed for Topology B, there is the additional root 16ũ + (4 +ṽ) 2 , which in the variables (A.2) is u(u + 8uv + 16(1 + u)v 2 ). We note that this new root appears together with all the preceding square roots and it is not possible to split the set of coincident square roots in sub-sets that contain a smaller number of roots at the time. Therefore we are forced to linearize 16ũ + (4 +ṽ) 2 together with all the preceding four roots. However, the diophantine equation associated with the fifth root 16ũ + (4 +ṽ) 2 cannot be solved. The reason is that it is associated with a polynomial of degree 8 in the variable z and degree 20 in the variable w, and it was not possible to find a parametrization, such as (A.18), which lower the degree of the polynomial equation so that it admits rational solutions.
Our alphabet in the new variables is constituted by linear letters and a single square root that remains in some of the repeated integrations.