A non-planar two-loop three-point function beyond multiple polylogarithms

We consider the analytic calculation of a two-loop non-planar three-point function which contributes to the two-loop amplitudes for $t \bar{t}$ production and $\gamma \gamma$ production in gluon fusion through a massive top-quark loop. All subtopology integrals can be written in terms of multiple polylogarithms over an irrational alphabet and we employ a new method for the integration of the differential equations which does not rely on the rationalization of the latter. The top topology integrals, instead, in spite of the absence of a massive three-particle cut, cannot be evaluated in terms of multiple polylogarithms and require the introduction of integrals over complete elliptic integrals and polylogarithms. We provide one-fold integral representations for the solutions and continue them analytically to all relevant regions of the phase space in terms of real function, extracting all imaginary parts explicitly. The numerical evaluation of our expressions becomes straightforward in this way.


Introduction
The study of the mathematical properties of multiloop Feynman integrals has received increasing attention in the last years both by the physics and the mathematics communities. The high precision reached by the experimental measurements carried out at the LHC, in fact, requires on the theory side the calculation of multi-(typically two-or three-) loop corrections to various complicated 2 → 1, 2 → 2 and 2 → 3 processes, including the exchange of massless and massive virtual particles. Indeed, unravelling their mathematical structures will be crucial to handle the complexity of such calculations and, on the long run, it could help gain a deeper understanding of quantum field theory itself. Impressive progress has been achieved in the last years in both directions. The development of new techniques for the calculation of multiloop Feynman integrals has rendered many previously out-ofreach calculations feasible. Among these, a fundamental role has been played by integration-by-parts reductions [1][2][3] and the differential equation method [4][5][6], augmented more recently by the use of a canonical basis [7,8]. A posteriori, a reason for the success of these techniques can be traced back to the almost concurrent "discovery" that the ǫ expansion of a rather large class of multiloop (mainly massless) Feynman integrals can be expressed in terms of a very well understood class of special functions, dubbed multiple polylogarithms [9][10][11][12]. In particular the understanding of their algebraic and analytical properties, also thanks to the development of the symbol and coproduct formalism [13][14][15][16][17], together with stable routines for their evaluation on the whole complex plane [18][19][20], have played a crucial role in the last years. Nevertheless, as it has been known for a long time, starting at the two-loop order the ǫ expansion of Feynman integrals can involve new mathematical structures which lie beyond the realm of multiple polylogarithms. The first and best studied example is that of the two-loop massive Sunrise graph [21][22][23][24][25][26][27][28][29][30][31][32][33][34]. Its imaginary part in d = 4 space-time dimensions is known to be expressible as linear combination of complete elliptic integrals of first and second kind. The real part of the graph can be reconstructed by a dispersive relation and involves integrals over elliptic integrals. In the last years quite some effort has been spent on studying the new functions appearing in the Sunrise graph, trying to extend the properties of multiple polylogarithms to embrace also elliptic generalizations of the latter [27,[29][30][31][32][33]. Much progress has been made but a complete understanding of these functions is still missing, in particular regarding their analytic continuation and numerical evaluation. Moreover, it is not clear whether they can be easily used to express also graphs with more complicated kinematics (i.e. three-and four-point functions). Last but not least, the interplay of these functions with the simpler multiple polylogarithms and their iterative structure remains unclear in context with the method of differential equations. There, simpler integrals appear as inhomogeneous terms in the differential equations of more complicated integrals. The solution of these equations requires on the one hand the solution of the homogeneous equations and, on the other, the integration over the inhomogeneous piece. Multiple polylogarithms (and in general Chen iterated integrals [35]) are particularly well suited for this scope, as they can be defined as iterative integrals over a set of differential forms. Of course, it is a priori unclear whether such a feature can be expected to hold for arbitrarily complicated Feynman graphs.
More recently, a new picture has started to emerged, where a generalization of this approach to more complicated cases becomes possible. When dealing with a system of coupled differential equations, the most non-trivial step consists in solving the homogeneous part of the system. In [36] it was shown that, for arbitrarily complicated cases, integral representations of the homogeneous solutions can be found by computing the maximal cut of the graph under consideration. Once the homogeneous solution is known, an integral representation for the inhomogeneous solution is provided by Euler's variation of constants. Even if usually such integrals cannot be expressed in terms of known special functions, from a practical point of view we are particularly interested in obtaining results that allow for fast and reliable numerical evaluations over the physical phase space. In this context, it was shown in [34] that the study of the imaginary part and of the corresponding dispersion relations of Feynman integrals within the differential equations framework can facilitate obtaining compact one-fold integral representations for the two-loop massive Sunrise and the Kite integral 1 . There are indications that more complicated Feynman integrals with different and unrelated kinematics can be casted in a similar form, see for example [38,39] and more recently a set of two-loop planar double boxes relevant for H+jet production [40].
The new ideas summarized above provide us with the tools required to start successfully studying more examples of relevant Feynman diagrams whose ǫ expansions do not evaluate to multiple polylogarithms. In this paper, we focus on a three-point non-planar topology which is relevant for the two-loop QCD corrections to the processes gg → tt and gg → γγ through a massive top loop. Similar integrals appear in the non-planar sectors of H+jet and HH production. The rest of the paper is organized as follows. In Section 2 we describe the problem and establish the notation. We continue in Section 3 showing how the ǫ expansion of the simpler subtopologies can be integrated in terms of multiple polylogarithm with a new algorithm based on the differential equations in canonical form. In Section 4 we consider the top topology, which cannot be integrated in terms of multiple polylogarithms. We show how to solve the coupled system of differential equations satisfied by its two master integrals based on information extracted from the maximal cut; we then compare it to the result obtained by solving directly the second order differential equation. We thoroughly study the required functions in Section 4 and Appendix B, where we also show how to analytically continue our solution to all relevant regions of the phase space. Finally in Section 4.3 we use the results of the previous sections to build up the inhomogeneous solution for the top topology in analytic form. We write the results as one-fold integrals over rational and irrational functions, multiple polylogarithms and complete elliptic integrals of first and second kind. We then draw our conclusions in Section 5.

Differential equations
We consider the Feynman integrals family I a1,a2,a3,a4,a5,a6,a7 defined as q p 1 p 2 = I a1,a2,a3,a4,a5,a6,a7 , and ǫ = (4 − d)/2. We stress that the triangle topology above is obtained by limiting ourselves to negative powers of the last propagator, i.e. a 7 < 0. Four of the six propagators are massive (thick lines), while only one of the three external legs is off-shell, i.e. p 2 1 = p 2 2 = 0 and q 2 = (p 1 + p 2 ) 2 = s . Integrals in family (2.1) with different values of the exponents a i may be reduced to master integrals using integration by parts reductions as implemented for example in Reduze 2 [41][42][43][44]. As a result we find 11 different master integrals: 9 for the subtopologies and 2 for the top topology.
As we will see explicitly, the nine master integrals of the subtopologies can be expressed in terms of multiple polylogarithms only. For them we choose a canonical basis using the method described in [45]. 2 The remaining two integrals, instead, contain elliptic generalizations of the multiple polylogarithms and, therefore, it is not possible to find for them a canonical basis as defined in [8]. 3 This can be motivated using the ideas presented in [49]. There it was shown for various examples that the possibility of decoupling the differential equations in the limit d → 4 (and therefore of writing the result in terms of multiple polylogarithms) is signaled by a degeneracy of the integration by parts identities in the limit of even numbers of dimensions, d → 2 n with n ∈ N.
We first consider the part of the Euclidean region where s < −4 m 2 and employ the following basis  where we added the number of different denominators and Reduze's sector id as a superscript to the I integrals for easier identification.
In order to write down the system of differential equations we introduce the massless ratio We introduce the vector m = (m 1 , . . . , m 9 ) for the master integrals of the subtopologies, which fulfil canonical differential equations. In matrix notation we obtain with the letters and the matrices The root appearing in the normalization of m 4 leads to l 2 and l 3 , while the root appearing in the normalization of m 6 leads to l 4 and l 5 . We note that, as long as we limit ourselves to m 1 , . . . , m 9 , these pairs of roots never mix and we could rationalize them separately with two different changes of variables. The two structures, nevertheless, mix up once considering the differential equations for the two master integrals of the top sector, m 10 and m 11 . With the basis given in (2.3), the differential equations for the master integrals of the top topology read d dx where B(x) and D(x) are two 2 × 2 matrices that do not depend on ǫ, while N 10 (ǫ; x) and N 11 (ǫ; x) contain the dependence on the simpler subtopologies. The matrices of the homogeneous part are while the non-homogeneous terms read From the latter it is clear that in the top topology all letters appear simultaneously. In our normalization the integrals m j = m j (ǫ; x) (j = 1, . . . , 11) have a Taylor series expansion at ǫ = 0, and we wish to calculate the expansion coefficients up to n = 4.

Integration of the subtopologies
We first consider the evaluation of the subtopologies m 1 , . . . , m 9 . Inserting the expansion (2.10) into the differential equation (2.4) and comparing coefficients of powers in ǫ gives which allows us to solve fully decoupled systems of differential equations order by order n bottom-up. The differential equations (3.1) have the form of the differential equations of multiple polylogarithms with root-valued letters. We find it useful to solve the full vector m in a uniform setup, where we consider all letters at the same time. This prevents a rationalization of the x dependence and therefore a traditional integration of the differential equations. Instead, we construct an ansatz for the solution using suitable target functions and require it to fulfil the differential equation (3.1). In other words, we integrate the symbol of the loop integrals. The basis of our function construction is the Duhr-Gangl-Rhodes algorithm [16], which we extend here to the case of an irrational alphabet. For simpler cases it is typically not strictly necessary to work with a canonical basis in order to be able to integrate the differential equations, provided they (partially) decouple. In contrast to this, our strategy here crucially relies on the differential equation (3.1) being in canonical form, such that we can read off the symbol letters of the solution. A similar approach for the integration of the differential equations has been used in [40] to calculate solutions through to weight two.
To construct our solutions in the region s < −4m 2 through to weight four we choose the functions with suitable arguments. We look for functions which do not introduce symbol letters beyond the alphabet (2.5) of the differential equation (3.1). For the arguments of ln we can just choose the original letters themselves. A classical polylogarithm Li j (z) has symbol entries z and 1 − z. Choosing z to be a power product of our original letters, c a0 l a1 1 · · · l a5 5 , with c a rational number, ensures z to factorize over our alphabet. Here, it is enough to consider a n ∈ and c a0 = ±1. We stress that this is too constraining in general. In particular, depending on the (choice of) letters it may be necessary to also consider roots of letters for factorization and allow the a n to be non-integer rational numbers. Note that there is an ambiguity related to the freedom of redefining the letters by expressions of the form (3.3). In our case, the definition of our letters absorbs any further explicit roots. In addition to the argument z of a polylogarithm we also require 1 − z to factorize in the sense (3.3) over our original alphabet. This condition is challenging to check with computer algebra systems, in particular for multivariate applications. We therefore resort to decompositions found by inserting numbers for the variables and applying heuristic integer relation searches to the logarithms of the involved factors [55][56][57] (more details will be given in [58]). Our heuristic search produces this admissible set: which is closed under z → 1 − z and z → 1/z, where z is an element of the set. For Li 2,2 (z 1 , z 2 ) we need to make sure that the symbol factors z 1 , 1 − z 1 , z 2 , 1 − z 2 and z 1 − z 2 do not introduce new letters. Our construction proceeds along the lines described above, but we do not list the generic result for the argument set here since it is too lengthy. We concentrate here on the part of the Euclidean region where s < −4m 2 and require all functions to be real valued for x > 4. Indeed, we find that all of these constraints can be satisfied and a solution in terms of such functions can be found, which satisfies the differential equations. Our functions read In order to fix the integration constants, we use the massive tadpole and the massless bubbles as input as well as various regularity conditions. We find very compact results for the final solutions of the subtopology integrals in terms of our choice of functions. These results will be a necessary building block to construct the full solution of the top topology. The explicit form of the results can be found in Appendix A and, in electronic format, in the ancillary file of the arXiv submission of this paper.

Homogenous solution
Having the expressions for all subtopologies in terms of multiple polylogarithms, we can now proceed and study the differential equations (2.7) for the top sector. We start by expanding in ǫ the inhomogeneous terms (2.9), up to ǫ 4 for j = 10, 11 .
10 (x) = 0 , N 10 (x) = 0 , N 10 (x) = 0 N where the letters l j have been defined above (2.5). We can now expand the right-and left-hand-side of (2.7) and we find, at any order n, Comparing with (4.2), we see that we get a non trivial non-homogeneous term only at order ǫ 4 , which is also the maximum order in the expansion we are interested in. In fact, the integrals I 6,63 1,1,1,1,1,1,0 and I 6,63 1,2,1,1,1,1,0 are finite, which can be checked e.g. with the methods of [59,60] implemented in Reduze 2.1, such that the first non-vanishing coefficients of m 10 and m 11 occur at order ǫ 4 .
The first step to solve (4.3) is to find a solution of its homogeneous part, which at every order in ǫ is a coupled system of two first-order linear differential equations d dx where we use a capital letter M instead of m to describe the homogeneous part of the solution and drop the superscript (n), since the form of the equation does not depend on n. In order to solve the system we need to find two independent set of solutions, i.e. a 2 × 2 matrix The inverse of the matrix G(x) reads where W (x) is the Wronskian of the solutions, is known, one can perform at every order in ǫ the rotation   where c (n) 10 and c (n) 11 are two integration constants and y 0 is a suitable boundary for the integration. In order to find an explicit solution for the matrix G(x) it is useful to recast (4.4) as a second order differential equation for the first master integral A general way to solve this equation is to compute the maximal cut of M 10 (x) as explained in detail in [36]. In that reference it was shown that the maximal cut of M 10 (x) computed on a specific integration contour can be written (up to an irrelevant overall numerical constant) as where we introduced the complete elliptic integral of the first kind defined as , for z ∈ C and ℜ(z) < 1 . It is straightforward to verify that (4.12) satisfies (4.11). A second independent solution can be obtained integrating on an independent contour. As explained in [36], this can be avoided in the case of elliptic integrals, as the second solution is simply given by (4.14) The two solutions (4.12) and (4.14) completely determine the first row of the matrix (4.5). Before proceeding, it is instructive to try to solve Eq. (4.11) directly, since in this case the equation is particularly simple. We define the auxiliary function which fulfills the differential equation Equation (4.16) is now in standard form and its solution can be easily written in terms of the complete elliptic integral of the first kind (4.13). Since the differential equation (4.16) has regular singular points in x = 0, x = 16 and x = ±∞, we must consider the solution of the latter in the three intervals 0 < x < 16, 16 < x < ∞ and −∞ < x < 0 . We start by considering the region 0 < x < 16 i.e. −16 m 2 < s < 0. This region is non physical and the master integrals m 10 , m 11 are real there. Following the same procedure, we can build up similar solutions in the remaining regions, i.e. the part of the Euclidean region where −∞ < s < −16 m 2 and the physical region s > 0. We derive explicitly those solutions in Appendix B. We use the notation I for the two pairs of solutions valid in the region a < x < b, and consequently G (a,b) (x) for the matrix of the solutions in the same region. The general solution of the homogeneous second order differential equation for 0 < x < 16 can be written as which determines the first row of the matrix (4.5) as Eqs. (4.18) are apparently different from the solutions determined by studying the maximal cut (4.12) and (4.14). Indeed, as the four functions are all solutions of the same second order homogeneous differential equation, a relation among them must exist such that only two functions are really independent. Indeed, it is well known that the elliptic integral of the first kind K(z) and its complementary K(1 − z) fulfil the following relations where we use the prescription z → z + i 0 + for definiteness. 4 Applying Eq. (4.19) twice to the solutions (4.12) and (4.14), we construct two more equivalent sets of solutions It is simple to check by direct insertion that these functions do indeed solve our second order differential equation. This proves the equivalence of the solutions found through the maximal cut with the ones found by solving directly the differential equation. The existence of more than one representation for the same solutions in terms of elliptic integrals of different arguments will turn out to be very important to write down the analytic continuation of the solution from 0 < x < 16 to the whole phase space in terms of explicitly real solutions, as explained in Appendix B.
We can now proceed with the solution of the system (4.3) in the region 0 < x < 16; here the most convenient representation is provided by the choice (4.18), which we will adopt from now on. In order to determine the second row of the matrix G(x) we plug (4.18) into (4.6) and find for consistency that where we introduced the complete elliptic integral of the second kind We write therefore the matrix of the solutions valid for 0 < x < 16 as . (4.24) We can now compute the Wronskian of the solutions. We find which is a direct consequence of the Legendre identity among the first two complete elliptic integrals (4.26) Alternatively, we can determine W (x) also without resorting to the Legendre identity. The fact that the Wronskian must be a linear function of x can be easily seen taking its derivative as with c an arbitrary integration constant. 5 The value of the constant c can be then fixed by computing the Wronskian for a fixed value of x. For example, we can study the behaviour of the solutions at the two boundaries, i.e. x → 0 + and x → 16 − . We find respectively, keeping only the leading behaviour, (2)) , By using (4.29) and (4.30) it is easy to verify explicitly that which fixes the constant c.

Analytic continuation of the inhomogeneous term
With the results of the previous sections we can solve the inhomogeneous differential equations for the top topology using Euler's variation of constants. In order to employ explicitly real building blocks, we need to analytically continue the homogeneous solutions and the inhomogeneous term to the various regions of the phase space. Details for the continuation of the homogeneous solutions are worked out in Appendix B. For the analytic continuation of the inhomogeneous term we give the explicit results in this section to highlight the different emerging functions. For clarity, in every region a < x < b we extract explicitly the imaginary part (whenever there is one) and introduce the functions R (a,b) (x) and Q (a,b) (x) respectively for the real and imaginary parts of N (4) For the analytical continuation to other kinematical regions we assign an infinitesimal positive imaginary part to s which translates to a negative imaginary part for x, x → x − i 0 + . We will also make use of letters from the alphabet depending on the region of phase space. For the region 4 < x < 16 we obtain R (4,16) (x) = 5 ln 2 (l 2 ) − 3 l 1 ζ 2 /2 + ln 2 (l 4 ) + Li 2 (−1/l 2 4 ) l 5 , We remark that x = 16 is actually no special point for the subtopologies and a single analytic expression is sufficient for the entire range x > 4. In the region 0 < x < 4, N 11 (x) stays real and we obtain R (0,4) (x) = 5 ln 2 (l 2 ) + 3 l 1 Cl 2 (−2 arccsc(2/l 1 )) l ′

5
, where the Clausen function Cl 2 (θ) = − θ 0 ln 2 sin t 2 dt arises from the dilogarithm of a pure phase factor, Cl 2 (θ) = ℑ Li 2 (e iθ ) , which can be seen from l 2 . In the region −4 < x < 0, the inhomogeneous terms develops an imaginary part and can be conveniently expressed in terms of where the first term of the real part arises from a purely imaginary ln(l 2 ) since l 2 2 = ( is a pure phase factor. Finally, for −∞ < x < −4 we obtain for the real and the imaginary part, respectively.

The inhomogeneous solution
We are now finally in the position of writing down the complete solution of (2.7) in all relevant regions of the phase space. For definiteness, we start in the Euclidean region 0 < x < 4, −4m 2 < s < 0, and then use the results of the previous section and of Appendix B to continue the solution to the other regions. In Appendix B, in particular, we build up matrices of solutions whose entries are explicitly real for a < x < b, and we show explicitly how to match the homogeneous solution moving from one region to the other. We recall that the expansions of m 10 and m 11 start at order ǫ 4 . We first write down the solutions form where (a, b) is either (0, 4) or (4,16), depending on whether y < 4 or y > 4. The integration constants must be fixed with two proper boundary conditions. Imposing that the integrals be regular as x → 16 − and go to zero for x → 0 + , we find c 10 = c 11 = 0  Note that now all integrals are explicitly real and all imaginary parts are explicit.
We can then go further and continue beyond the threshold for the production of two massive particles, i.e. s > 4m 2 , x < −4. Using again formulas (4.36), (B.16) and (B.17) we obtain m (4)  They can be evaluated with high precision starting directly from this integral representation.

Numerical evaluation
We verified our formulas in the previous section with a thorough comparison against SecDec 3, both in the physical and in the non-physical region. Numerical results for the top-level master integrals are shown in Figure 1. We would like to point out that our representation allows for particularly fast and precise evaluations also in the physical region of phase space. As an example for a result with higher precision, we obtain m 4 I 1,1,1,1,1,1,0 s=5m 2 ,d=4 ≈ −0.07776462028160023644086669458011467822536257409024 for one of our six-line master integrals.

Conclusions
Feynman integrals which evaluate to classes of functions outside the realm of multiple polylogarithms constitute the bottleneck for many multiloop calculations relevant for LHC phenomenology. To this day, only a very limited number of examples have been considered in the literature. An important step towards a more complete understanding of the new mathematical structures therefore consists of the study of more explicit examples in order to expose their common properties. In this paper we have considered the calculation of a two-loop non-planar three-point function with four massive internal propagators and one external off-shell leg using the method of differential equations. The differential equations for the subtopologies can be put in canonical form and integrated in terms of multiple polylogarithms over an irrational alphabet. In general, in the presence of different roots, traditional approaches to integrate the result in terms of Goncharov's polylogarithms fail, unless one is able to find a change of variables which linearizes or at least rationalizes the letters of the alphabet. We employed a new algorithm to perform the integration without any rationalization of the alphabet. This allowed us to write the results up to weight four in terms of simple logarithms, classical polylogarithms and Li 2,2 functions only. We moved then to the two master integrals of the top level sector, which fulfil two coupled differential equations in the momentum transfer squared. By studying the maximal cut of the two master integrals we showed how to solve the differential equations in terms of complete elliptic integrals of first and second kind; we then used Euler's variation of constants in order to write the solution as a one-fold integral over elliptic integrals and multiple polylogarithms. Finally, we used the same techniques described in [34] in order to continue our results to the Minkowski physical region. In this way we provided expressions which can be evaluated in a fast and reliable manner in the entire physical region of the phase space. We compared our results with SecDec 3 and found perfect agreement. We expect that similar techniques can be extended to study also more complicated three-and four-point functions, as in part confirmed by the results recently obtained in [40].

Acknowledgments
We would like to thank Stefan Groote, Ettore Remiddi, Robert Schabinger and Erich Weihs for valuable discussions. We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support.

A Solutions for the subtopologies
In this Appendix we report the explicit results for the subtopologies in the Euclidean region with s < −4m 2 . Having a canonical form, we can write the ǫ expansion through to weight four in a relatively compact form in terms of simple functions, including only logarithms, classical polylogarithms and Li 2,2 functions with arguments built from the alphabet defined in (2.5). Expanding the master integrals according to As discussed above, the Wronskian of the solutions must be a linear function of x, see (4.28). It is simple starting from (B.6) and (B.7) to prove that we have also for our solutions in this region.
Euclidean region 16 < x < ∞ The same idea can be applied in the remaining part of the Euclidean region. Similar to the Minkowski region, we start here with solutions constructed as complex linear combinations of the original functions which evaluate to real numbers for 16 < x < ∞,