Two-loop Master Integrals with the Simplified Differential Equations approach

We calculate the complete set of two-loop Master Integrals with two off mass-shell legs with massless internal propagators, that contribute to amplitudes of diboson $V_1V_2$ production at the LHC. This is done with the Simplified Differential Equations approach to Master Integrals, which was recently proposed by one of the authors.


Introduction
At the first run of the LHC, proton collisions have reached the record-setting high energies of 8 TeV and delivered luminosity of ∼ 23.3 fb −1 . At the second run of the LHC, which is expected to start in 2015, the energy and luminosity will be pushed much higher. In order to keep up with the increasing experimental accuracy as more data is collected, more precise theoretical predictions and higher loop calculations will be required.
With the better understanding of reduction of one-loop amplitudes to a set of Master Integrals (MI) based on unitarity methods [1,2] and at the integrand level via the OPP method [3,4], one-loop calculations have been fully automated in many numerical tools (some reviews on the topic are [5,6]). In the recent years, a lot of progress has been made towards the extension of these reduction methods to the two-loop order at the integral [7,8] as well as the integrand [9][10][11][12] level. Contrary to the MI at one-loop, which have been known for a long time already [13], a complete library of MI at two-loops is still missing. At the moment this seems to be the main obstacle to obtain a fully automated NNLO calculation framework similar to the one-loop case, that will satisfy the anticipated precision requirements at the LHC [14].
Starting from the works of [15][16][17], there has been a building consensus that the socalled Goncharov Polylogarithms (GPs) form a functional basis for many MI. A very fruitful method for calculating MI and expressing them in terms of GPs is the differential equations (DE) approach [18][19][20][21][22], which has been used in the past two decades to calculate various MI at two-loops [23][24][25][26][27][28]. In [29] a variant of the traditional DE approach to MI was presented, which was coined the Simplified Differential Equations (SDE) approach. In this paper a review of the SDE approach and an application to the four-point two-loop planar and non-planar MI with two different off-shell external legs and massless internal propagators is presented.
Those double-box integrals are needed in particular in order to compute NNLO QCD corrections to pp → V V * processes at LHC, where V V * are massive electroweak vector bosons. In particular the case of two different masses can be helpful to consider ZZ * or W + W − off-shell production as well as different bosons production like W Z (or V γ * , ZH, γ * γ * , etc.). Another significant improvement is related to the fact that pp → V V * is the main contribution to the background for testing Higgs boson decay into two vector bosons. A better understanding of this background could lead to improved measurements of Higgs couplings to vector bosons as well as the Higgs decay width [30,31].
For these reasons, recently two loop MI with two off-shell external legs have been intensively studied and calculated by different groups [27,28,32,33]. Furthermore those results have been already implemented to provide NNLO predictions for ZZ production at the LHC [34]. In this paper we calculate the same MI as those given in [27,28], namely the full set of two loop MI with different masses which contains both planar as well as non-planar Feynman diagrams, by using the independent SDE approach presented in [29].
The paper is organized as follows. In the Section 2 we explain and review the SDE method. In Section 3 we give the parameterization and notation of the variables describing the two-loop MI. As we will show, the same parameterization of the external legs can be used to parametrize both the planar and non-planar MI. In Section 4 we discuss our results, we show how these results can be analytically continued in order to accommodate both the Euclidean as well as the physical kinematical regions and compare with results obtained by direct numerical integration as well as with the analytic calculations given in [27]. We conclude in Section 5 and provide an overview of the topic and some perspective for future developments. In the Appendix we give few characteristic examples on how the boundary conditions are properly reproduced in our approach by the DE. Finally in the attached ancillary files, we provide our analytic results for all two-loop MI in terms of Goncharov polylogarithms together with explicit numerical results.

Simplified differential equations
The MI in this paper will be calculated with the SDE approach [29] and therefore in this Section a brief review is given of how the method is used in practice. Assume that one is interested in calculating an l−loop Feynman integral whose graph with external momenta {p j }, considered incoming, is shown on the left hand side in Figure 1. Throughout this paper the internal propagators are assumed massless. The Feynman integral on the left hand side in Figure 1 is a subset of the following class of l−loop integrals: where the denominators are defined in such a way that all scalar product invariants can be written as a linear combination of them. The exponents a i are integers and may be negative in order to accommodate non-trivial numerators. Any integral G a 1 ···an may be written as a linear combination of the MI G M I ({s ij }, ) with coefficients depending on the independent scalar products, s ij = p i · p j , and space-time dimension d, by the use of integration by parts (IBP) identities [35,36]. In the traditional DE method, the vector of MI G M I ({s ij }, ) is differentiated with respect to p i · ∂ ∂p j , and the resulting integrals are reduced by IBP to give a linear system of DE for G M I ({s ij }, ) [18,21]. The invariants, s ij , are then parametrised in terms of dimensionless variables, defined on a case by case basis, so that the resulting DE can be solved in terms of GPs. Usually boundary terms corresponding to the appropriate limits of the chosen parameters have to be calculated using for instance expansion by regions techniques [37,38].
In the SDE approach [29] a dimensionless parameter x is introduced so that the vector of MI G M I ({s ij }, ; x), which now depends on x, satisfies a system of differential equations in one independent variable. The expected benefit is that the integration of the DE naturally captures the expressibility of MI in terms of GPs and more importantly makes the problem independent on the number of kinematical scales (independent invariants) involved. As shown on the right hand side in Figure 1, the external incoming momenta are now parametrized linearly in terms of x as p i (x) = p i + (1 − x)q i , where the q i 's are a linear combination of the momenta {p i } such that i q i = 0. If p 2 i = 0, the parameter x captures the off-shellness of the external legs. The class of Feynman integrals in (2.1) are now dependent on x through the external momenta: Note that as x → 1, the original configuration of the loop integrals (2.1) is reproduced, which correspond to a simpler one with one scale less. We briefly summarize how the DE (2.2) may be solved in practice. The MI with least amount of denominators m 0 (m 0 = 3 in the case of two loops) are typically two-point integrals that satisfy homogeneous equations and can be easily calculated analytically by other methods. Furthermore, the form of 2.2 is such that MI with m denominators only depend on MI with at most m denominators. This structure of the DE makes it possible to first solve the MI with m 0 + 1 denominators, then those with m 0 + 2 denominators and so forth. In other words, in practice the DE may be solved in a bottom-up approach.
Assume that all MI with m ≤ m denominators are known and already expressed in the desired form, the meaning of which will become clear below. The DE of MI with m + 1 denominators can be written schematically in the form: The equation could now be straightforwardly expanded in and integrated, provided that the right-hand side is free of singularities at x → 0. In order to deal with possible singularities at x = 0 we need therefore to determine the re-summed part of the solution in this limit.
The righthand side of (2.5) can be rewritten schematically as a sum of a singular and regular term at x = 0: The singular term is integrated exactly and the solution of (2.5) becomes: where the first term C({s ij }, ) is a constant in x but may be dependent on the the kinematical invariants. The rightmost term in Eq.(2.7), is safely expanded in and expressed in terms of GPs. To this end the integrating factors M in Eq.(2.5) should be rational functions of x in the limit → 0. This is a sufficient condition for the chosen x-parametrisation to result in a differential equation solvable in terms of GPs.
As was noted in [29], in the SDE approach the boundary terms, in all cases studied until now including one and two loops MI, are always naturally captured by the integrated singularities in the SDE (2.7) themselves at x = 0, which is precisely the lower integration boundary of the GPs. In other words, in all cases studied the integrated singular terms in (2.7) correctly describe the behaviour of M G m+1 as x → 0 and thus the constant C({s ij }, ) vanishes. In this way, the SDE method is well suited for directly and efficiently expressing the MI in terms of GPs without the need for an independent evaluation of the MI at the boundary x = 0. In fact all re-summed parts in the limit x → 0 are fully determined by the one-scale MI involved in the system -in the case of the present calculation, two-point integrals, two three-point integrals and double one-loop integrals -that satisfy homogenous differential equations. These integrals are simple enough to be evaluated by other methods and in fact they are all known for the cases of interest. As we have verified using the expansion by regions technique this property is closely related to the x-parametrization that allows the expression of the limit x → 0 of all integrals involved in terms of the one-scale integrals, see Appendix A.
The integrand of the remaining integral in (2.7) is regular at the boundary x = 0. After expanding in and performing partial fraction decomposition in x, the integral is directly expressible in terms of Goncharov polylogarithms, which are an iterative class of functions that generalize the usual logarithms and polylogarithms [15,17]: where in general α i , x ∈ C. Note that the integration boundary in (2.7) was chosen to be x = 0 precisely in order to directly express the integrals in terms of GPs. If another boundary point x = x 0 is chosen, the integrals in (2.7) result in slightly more complicated expressions made up of differences of GPs. We will denote the first set of arguments of the GPs, i.e. those contributing to the weight, as indices. The SDE method therefore directly expresses the MI in terms of a well defined functional basis of GPs, independently on the number of kinematical scales involved. When the DE are coupled, the homogeneous factor H is a matrix in (2.4). For all cases considered so far, proceeding the same way as in the uncoupled case, the diagonal elements are used to determine the integrating factors and the singularities at x = 0 in the inhomogeneous terms to determine the re-summed part of the solutions. The homogeneous matrix of the reduced system of DE is then a strictly triangular matrix at order 0 and the system becomes effectively uncoupled, and may be easily integrated order by order in as was for example explained in [29]. One may also envisage to solve the reduced system by using a Magnus or Dyson series expansion [39,40]. The latter approach is possible because the non-diagonal piece of the 0 matrix is stricly triangular and thus nilpotent. We refer to [33] for a similar discussion for the case of coupled DE.
In few cases singularities of the form x −2+l appears in the inhomogeneous terms. In these cases, by using the transformation x → 1/x we are able to bring the DE in their normal form with singularities of the form x −1+l and then follow the same procedure explained Figure 3. The parametrization of external momenta for the three planar double boxes of the families P 12 (left), P 13 (middle) and P 23 (right) contributing to pair production at the LHC. All external momenta are incoming. above. In fact as can be seen in the Appendix A, the DE after the transformation do correctly reproduce the boundary conditions in terms of one-scale integrals as usual. Of course the results are expressed back in terms of x by applying the x → 1/x transformation on the argument of the GPs. As far as the parametrization of the external momenta in x is concerned, we noticed that for all the cases that we considered it was enough for us to choose the parametrization of the external legs such that after pinching internal lines the resulting triangles with three off-shell legs, if they appear, have the form given in Figure 2 [29].

Notation and x-parametrization of the Master Integrals
We are interested in calculating the MI of two-loop QCD amplitude corrections contributing to diboson pair production at the LHC, where the two outgoing particles may have different masses: The colliding partons have massless momenta q 1 , q 2 and the outgoing particles have massive momenta q 3 , q 4 . Both the planar and non-planar diagrams contributing to (3.1) have already been calculated with the traditional DE method [27,28]. There are in total six families of MI whose members with the maximum amount of denominators are graphically shown in Figures 3 and 4. Three of these, namely those shown in Figure 3, contain only planar MI and will therefore be referred to as the planar families in this paper. They will be denoted by P 12 , P 13 and P 23 as was done in [27] and contain 31, 29 and 28 MI respectively. The other three families, shown in Figure 4, contain planar MI with up to 6 denominators as well as non-planar MI with 6 and 7 denominators and will be referred to as the non-planar families of MI. These non-planar families will be denoted by N 12 , N 13 and N 34 as was done in [28] and contain 35, 43 and 51 MI respectively. In this paper we have used both FIRE [41] and Reduze 2 [42] to perform the IBP reduction to the MI. In order to parametrize the external momenta in terms of x as was dicussed in the previous Section we kept in mind Figure 2. In particular, the parametrization of P 12 and P 13 were chosen such as to satisfy the requirement of having off-shell triangles of the form shown in Figure 2 after pinching the internal line(s) between the two massless legs in Figures 3. The parametrization of P 23 was found by permuting the external legs accordingly 1 . All non-planar diagrams shown in Figures 4 may be transformed, after pinching enough internal lines, to off-shell triangles with the exact same external momenta as those off-shell triangles found after pinching the families P 12 and P 13 . Therefore in the SDE approach it is sufficient to choose the same parameterization of the external legs for the non-planar families as was chosen for the planar ones in order for the families to contain off-shell triangles of the form shown in Figure 2. The incoming external momenta for all families are therefore parametrized as follows: where the notation p i···j = p i +· · ·+p j is used. As x → 1, the external momentum q 3 becomes massless, such that (3.2) also captures the case where only one outgoing particle is offshell. The MI depend in total on 4 variables, namely the traditional Mandelstam variables S = (q 1 + q 2 ) 2 , T = (q 1 − q 3 ) 2 and (squared) particle masses M 2 3 = q 2 3 , M 2 4 = q 2 4 . The x-parameterization (3.2) of the external momenta results in these variables being related to the parameter x, the external off-shell mass squared q and the two independent scalar products s 12 and s 23 of our choice that were defined in (3.2): The other Mandelstam variable U = (q 1 − q 4 ) 2 is related to the other variables as usual via S + T + U = M 2 3 + M 2 4 . After finding the x-parameterization (3.2), the class of loop integrals describing the planar families P 12 , P 13 and P 23 are now explicitly expressed in x as: , (3.5) 1 Note that for the family P23 it is not possible to get off-shell triangles by pinching the internal lines.
where γ E is the usual Euler-Mascheroni constant. Similarly, the class of integrals for the non-planar families are: Using the notation given in Eq. (3.4), (3.5) and (3.6) for the class of planar loop integrals, the indices a 1 · · · a 9 for the list of MI in the planar families are as follows 2 : Most of the above planar MI also appear in the non-planar families, whose list of MI have the following indices a 1 · · · a 9 in terms of the notation of Eq.   , T < 0, U < 0, where q 2 ⊥ is the transverse momentum squared of each massive particle in the centre of mass frame. In terms of our variables the physical region is therefore: The second inequality in (3.17) causes a branching of the s 23 variable in the physical region: x > 1, s 23 < 0, s 12 + s 23 > q, q > 0 s 23 > 0, s 12 + s 23 < q, s 12 > q/x.

(3.18)
We note that the boundary terms of our MI are evaluated in the Euclidean region x = 0 as explained in the previous Section. In the next Section we will discuss our results for the MI and describe how they can be analytically continued from the Euclidean region x < 1 to the physical region where x > 1 (3.17).

Results for Master Integrals
With the chosen x-parametrization, the solutions of the DE are all expressed in terms of GPs (2.8). The set of all indices I(P i ) of the GPs appearing in the MI of each planar family is: Note that the parameter x is not part of the set of indices but appears in the argument of the GPs as we integrate the DE. This can be for example explicitly seen in the solution of the scalar double box of the family P 13 : The above MI is expressed in terms of a common factor A 3 ( ) As one notices from the example G P 13 011111011 above, the solutions are in general not uniform in the weight of the GPs [22]. The Goncharov polylogarithms can be numerically evaluated with the Ginac library [43]. We have not tried to simplify our analytical expressions by the use of symbol and co-product techniques [44,45], as we were mostly interested in showing the applicability of the SDE method, but this is expected to be straightforward. The full expressions for all MI are available in the attached ancillary files.
In order to compute the MI in arbitrary kinematics, especially in the physical region, the GPs have to be properly analytically continued. In general all variables, including the momenta invariants s ij and the parameter x would acquire an infinitesimal imaginary part, s ij → s ij + iδ s ij η, x → x + iδ x η, with η → 0. The parameters δ s ij and δ x are determined as follows: in the first place the input data to the DE, namely the two-point two-loop master integrals need to be properly defined in each kinematical region. For instance writing the two-point function G P 12 001000011 as implies constraints on the imaginary parts of the quantities involved. Moreover the second graph polynomial [46], F, in the limit η → 0, should acquire a (definite-)negative imaginary part. Combining these two classes of constraints on the infinitesimal parameters δ s ij and δ x , fixes the imaginary part of all the GPs involved and we have checked that the result for the MI is identical in the limit η → 0 and moreover agrees with the one obtained by other calculations.
We have performed several numerical checks of all our calculations. The numerical results have been compared with those provided by the numerical code SecDec [47][48][49][50] in the Euclidean region and with analytic results presented in [27,28] for the physical region. In all cases we find perfect agreement and reference numerical results can be found in the ancillary files.

Conclusions
In this paper we calculated the complete set of two-loop MI with massless propagators contributing to diboson production at the LHC and performed an independent calculation of the MI presented in [27]. The MI were calculated with the SDE method proposed in [29]. After choosing the x-parameterization as given in this paper, all MI are expressible in terms of GPs, independently on the number of kinematical scales involved, and moreover the boundary terms are directly captured by integrating the singularities at the boundary point x = 0. This seems to be the case for all MI calculated with the SDE method, which makes it an efficient and practical tool for the calculation of MI since it does not require the independent derivation of limits of those MI. We have shown that one x-parameterization of the external momenta is enough to parametrize both the planar as well as non-planar integrals and express them in terms of GPs with the argument x. The indices of the GPs are functions of the other three parameters that are linearly related to the scalar products of the external momenta when x → 1.
The phase space point at the boundary x = 0 lies in the Euclidean region and we discussed how to analytically continue the solutions for the MI to the physical region. All our results have been numerically compared with Secdec in the Euclidean region and with the analytical results of [27,28] in the physical region for various phase space points and perfect agreement has been found. All MI are attached as ancillary files to this paper together with reference numerical results for the MI with 7 denominators. Our results may also be used to calculate virtual two-loop QCD corrections contributing to processes at the LHC such as W Z, V γ * , ZH, or to calculate sub-diagrams of three-loop Feynman diagrams.
There are still several open issues which need further study. In particular it will be useful to find all x-parameterizations of the external momenta that a priori lead to DE directly expressible in terms of Goncharov polylogarithms. Moreover, as the DE are directly related to the IBP identities, it is a very intriguing question if such a parameterization can be determined by studying the IBP identities themselves. Another subject of future research is the extension of the SDE method to integrals with massive internal lines and more external legs.
(A.1) The integrating factor of the DE of I behaves as M (x) ∼ x 2 as x → 0 and we are therefore interested in the behaviour lim x→0 M I(x) = lim x→0 x 2 I = lim x→0 x 2 R∈R I (R) , where R denotes the set of all regions contributing. Performing the expansions for all the regions given above, it can be shown by power counting that we have a vanishing limit lim x→0 (x 2 I (R) ) = 0 for all regions except the purely soft region R = (S, S), which is defined by S = k ∼ xp 123 (k is either k 1 or k 2 ). Therefore, the behaviour of M I(x) as x → 0 is fully determined by the soft region: 2), as can be seen from the ancillary files. We conclude that the DE correctly captures the single-scale behaviour ∼ G P 13 001010001 as x → 0 via the re-summed term.
• Non-planar double box: G N 12

101111011
The integral equals: 3) The integrating factor of the DE of I behaves as M (x) ∼ x 4+2 as x → 0 and we are therefore interested in the behaviour lim x→0 M I(x) = lim x→0 x 4+2 R∈R I (R) . Performing the expansions for all the regions given above, it can be shown by power counting that lim x→0 (x 4+2 I (R) ) = 0 for all regions except the purely soft region R = (S, S). Therefore, the behaviour of M I(x) as x → 0 is fully determined by the soft region: The integral equals: (A.5) The integrating factor of the DE of I behaves as M (x) ∼ x 3 as x → 0 and we are therefore interested in the behaviour lim x→0 M I(x) = lim x→0 x 3 R∈R I (R) . In this case, the regions contribute to the behaviour of M I(x) as x → 0: