Maximal cuts and differential equations for Feynman integrals. An   application to the three-loop massive banana graph

Amedeo Primo; Lorenzo Tancredi

arXiv:1704.05465·hep-ph·November 26, 2018

Maximal cuts and differential equations for Feynman integrals. An application to the three-loop massive banana graph

Amedeo Primo, Lorenzo Tancredi

PDF

TL;DR

This paper develops a method using maximal cuts and elliptic integrals to solve the differential equations for three-loop massive banana Feynman integrals with equal masses, providing explicit solutions.

Contribution

It introduces a novel approach combining maximal cuts and elliptic integrals to explicitly solve the differential equations for the three-loop massive banana graph.

Findings

01

Homogeneous solutions expressed in terms of elliptic integrals.

02

All independent solutions obtained via contour integration of the maximal cut.

03

Inhomogeneous solutions derived using Euler's variation of constants.

Abstract

We consider the calculation of the master integrals of the three-loop massive banana graph. In the case of equal internal masses, the graph is reduced to three master integrals which satisfy an irreducible system of three coupled linear differential equations. The solution of the system requires finding a $3 \times 3$ matrix of homogeneous solutions. We show how the maximal cut can be used to determine all entries of this matrix in terms of products of elliptic integrals of first and second kind of suitable arguments. All independent solutions are found by performing the integration which defines the maximal cut on different contours. Once the homogeneous solution is known, the inhomogeneous solution can be obtained by use of Euler's variation of constants.

Equations441

\displaystyle\leavevmode\hbox{\raise-9.78711pt\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{\includegraphics[scale={0.8},trim=0.0pt 5.69046pt 0.0pt 0.0pt]{fig/fig2LSun.pdf} }}}}}}=S(d;u)=\int\,\frac{\mathfrak{D}^{d}k_{1}\,\mathfrak{D}^{d}k_{2}}{[k_{1}^{2}-m^{2}][k_{2}^{2}-m^{2}][(k_{1}-k_{2}-p)^{2}-m^{2}]}\,,

\displaystyle\leavevmode\hbox{\raise-9.78711pt\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{\includegraphics[scale={0.8},trim=0.0pt 5.69046pt 0.0pt 0.0pt]{fig/fig2LSun.pdf} }}}}}}=S(d;u)=\int\,\frac{\mathfrak{D}^{d}k_{1}\,\mathfrak{D}^{d}k_{2}}{[k_{1}^{2}-m^{2}][k_{2}^{2}-m^{2}][(k_{1}-k_{2}-p)^{2}-m^{2}]}\,,

[\frac{d ^{2}}{d u ^{2}} + (\frac{1}{u} + \frac{1}{u - 1} + \frac{1}{u - 9}) \frac{d}{d u} + (- \frac{1}{3 u} + \frac{1}{4 ( u - 1 )} + \frac{1}{12 ( u - 9 )})] S (u) = 0,

[\frac{d ^{2}}{d u ^{2}} + (\frac{1}{u} + \frac{1}{u - 1} + \frac{1}{u - 9}) \frac{d}{d u} + (- \frac{1}{3 u} + \frac{1}{4 ( u - 1 )} + \frac{1}{12 ( u - 9 )})] S (u) = 0,

Cut (S (u)) = \oint_{C} \frac{d b}{\pm b ( b - 4 ) ( b - ( u - 1 ) ^{2} ) ( b - ( u + 1 ) ^{2} )} = \oint_{C} \frac{d b}{\pm R _{4} ( b , u )},

Cut (S (u)) = \oint_{C} \frac{d b}{\pm b ( b - 4 ) ( b - ( u - 1 ) ^{2} ) ( b - ( u + 1 ) ^{2} )} = \oint_{C} \frac{d b}{\pm R _{4} ( b , u )},

0 < 4 < (u - 1)^{2} < (u + 1)^{2} .

0 < 4 < (u - 1)^{2} < (u + 1)^{2} .

\oint_{C_{\infty}} \frac{d b}{R _{4} ( b , u )} = 0,

\oint_{C_{\infty}} \frac{d b}{R _{4} ( b , u )} = 0,

\oint_{C_{\infty}} \frac{d b}{R _{4} ( b , u )} = \oint_{C_{1}} \frac{d b}{R _{4} ( b , u )} + \oint_{C_{3}} \frac{d b}{R _{4} ( b , u )} = 0,

\oint_{C_{\infty}} \frac{d b}{R _{4} ( b , u )} = \oint_{C_{1}} \frac{d b}{R _{4} ( b , u )} + \oint_{C_{3}} \frac{d b}{R _{4} ( b , u )} = 0,

\oint_{C_{1}} \frac{d b}{R _{4} ( b , u )} = 2 i

\oint_{C_{1}} \frac{d b}{R _{4} ( b , u )} = 2 i

\oint_{C_{3}} \frac{d b}{R _{4} ( b , u )} = - 2 i \int_{(u - 1)^{2}}^{(u + 1)^{2}} \frac{d b}{- R _{4} ( b , u )},

\oint_{C_{2}} \frac{d b}{- R _{4} ( b , u )} = 2 i (\int_{- \infty}^{0} + \int_{(u + 1)^{2}}^{+ \infty}) \frac{d b}{R _{4} ( b , u )} .

\oint_{C_{2}} \frac{d b}{- R _{4} ( b , u )} = 2 i (\int_{- \infty}^{0} + \int_{(u + 1)^{2}}^{+ \infty}) \frac{d b}{R _{4} ( b , u )} .

J (u) = \int_{0}^{4} \frac{d b}{- R _{4} ( b , u )} = \int_{(u - 1)^{2}}^{(u + 1)^{2}} \frac{d b}{- R _{4} ( b , u )},

J (u) = \int_{0}^{4} \frac{d b}{- R _{4} ( b , u )} = \int_{(u - 1)^{2}}^{(u + 1)^{2}} \frac{d b}{- R _{4} ( b , u )},

I (u) = \int_{4}^{(u - 1)^{2}} \frac{d b}{R _{4} ( b , u )} = (\int_{- \infty}^{0} + \int_{(u + 1)^{2}}^{+ \infty}) \frac{d b}{R _{4} ( b , u )} .

I_{C_{a}} = \oint_{C_{a}} \frac{d a}{( a - a _{1} ) ( a - a _{2} ) ( a - a _{3} ) ( a - a _{4} )} = \oint_{C_{a}} \frac{d a}{R ( a , a _{1} , a _{2} , a _{3} , a _{4} )},

I_{C_{a}} = \oint_{C_{a}} \frac{d a}{( a - a _{1} ) ( a - a _{2} ) ( a - a _{3} ) ( a - a _{4} )} = \oint_{C_{a}} \frac{d a}{R ( a , a _{1} , a _{2} , a _{3} , a _{4} )},

a_{1} < a_{2} < a_{3} < a_{4}

a_{1} < a_{2} < a_{3} < a_{4}

I (a_{1}, a_{2}) = \int_{a_{1}}^{a_{2}} \frac{d a}{- R ( a , a _{1} , a _{2} , a _{3} , a _{4} )} = \int_{a_{3}}^{a_{4}} \frac{d a}{- R ( a , a _{1} , a _{2} , a _{3} , a _{4} )},

I (a_{1}, a_{2}) = \int_{a_{1}}^{a_{2}} \frac{d a}{- R ( a , a _{1} , a _{2} , a _{3} , a _{4} )} = \int_{a_{3}}^{a_{4}} \frac{d a}{- R ( a , a _{1} , a _{2} , a _{3} , a _{4} )},

I (a_{2}, a_{3}) = \int_{a_{2}}^{a_{3}} \frac{d a}{R ( a , a _{1} , a _{2} , a _{3} , a _{4} )} = (\int_{- \infty}^{a_{1}} + \int_{a_{4}}^{+ \infty}) \frac{d a}{R ( a , a _{1} , a _{2} , a _{3} , a _{4} )} .

I (a_{1}, a_{2}) ⟶ t^{2} = \frac{( a _{4} - a _{2} ) ( a - a _{1} )}{( a _{2} - a _{1} ) ( a _{4} - a )},

I (a_{1}, a_{2}) ⟶ t^{2} = \frac{( a _{4} - a _{2} ) ( a - a _{1} )}{( a _{2} - a _{1} ) ( a _{4} - a )},

I (a_{2}, a_{3}) ⟶ t^{2} = \frac{( a _{1} - a _{3} ) ( a - a _{2} )}{( a _{3} - a _{2} ) ( a _{1} - a )},

I (a_{1}, a_{2}) = \frac{2}{( a _{3} - a _{1} ) ( a _{4} - a _{2} )} K (w_{1}),

I (a_{1}, a_{2}) = \frac{2}{( a _{3} - a _{1} ) ( a _{4} - a _{2} )} K (w_{1}),

I (a_{2}, a_{3}) = \frac{2}{( a _{3} - a _{1} ) ( a _{4} - a _{2} )} K (1 - w_{1}),

K (w) = \int_{0}^{1} \frac{d z}{( 1 - z ^{2} ) ( 1 - w z ^{2} )} \mbox w i t h ℜ (w) < 1,

K (w) = \int_{0}^{1} \frac{d z}{( 1 - z ^{2} ) ( 1 - w z ^{2} )} \mbox w i t h ℜ (w) < 1,

w_{1} = \frac{( a _{2} - a _{1} ) ( a _{4} - a _{3} )}{( a _{3} - a _{1} ) ( a _{4} - a _{2} )} .

w_{1} = \frac{( a _{2} - a _{1} ) ( a _{4} - a _{3} )}{( a _{3} - a _{1} ) ( a _{4} - a _{2} )} .

\displaystyle\leavevmode\hbox{\raise-20.98895pt\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{\includegraphics[scale={1.},trim=0.0pt 0.0pt 0.0pt -7.68236pt]{fig/fig3Lbubble.pdf} }}}}}}=I_{a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7},a_{8},a_{9}}\Big{|}_{a_{5},\cdots\,,a_{9}<0}

\displaystyle\leavevmode\hbox{\raise-20.98895pt\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{\includegraphics[scale={1.},trim=0.0pt 0.0pt 0.0pt -7.68236pt]{fig/fig3Lbubble.pdf} }}}}}}=I_{a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7},a_{8},a_{9}}\Big{|}_{a_{5},\cdots\,,a_{9}<0}

= \int \frac{D ^{d} k _{1} D ^{d} k _{2} D ^{d} k _{3} ( k _{3}^{2} ) ^{- a_{5}} ( k _{1} \cdot p ) ^{- a_{6}} ( k _{2} \cdot p ) ^{- a_{7}} ( k _{3} \cdot p ) ^{- a_{8}} ( k _{1} \cdot k _{2} ) ^{- a_{9}}}{[ k _{1}^{2} - m ^{2} ] ^{a_{1}} [ k _{2}^{2} - m ^{2} ] ^{a_{2}} [( k _{1} - k _{3} ) ^{2} - m ^{2} ] ^{a_{3}} [( k _{2} - k _{3} - p ) ^{2} - m ^{2} ] ^{a_{4}}},

\int D^{d} k = \frac{( m ^{2} ) ^{\frac{4 - d}{2}}}{i π ^{d /2} Γ ( \frac{4 - d}{2} )} \int \frac{d ^{d} k}{( 2 π ) ^{d}}

\int D^{d} k = \frac{( m ^{2} ) ^{\frac{4 - d}{2}}}{i π ^{d /2} Γ ( \frac{4 - d}{2} )} \int \frac{d ^{d} k}{( 2 π ) ^{d}}

\int \frac{D ^{d} k}{k ^{2} - m ^{2}} = \frac{2 m ^{2}}{d - 2} .

\int \frac{D ^{d} k}{k ^{2} - m ^{2}} = \frac{2 m ^{2}}{d - 2} .

I_{1} (ϵ; s) =

I_{1} (ϵ; s) =

I_{2} (ϵ; s) =

I_{3} (ϵ; s) =

I_{0} (ϵ; m^{2}) = I_{2, 2, 2, 0, 0, 0, 0, 0, 0} = 1 .

I_{0} (ϵ; m^{2}) = I_{2, 2, 2, 0, 0, 0, 0, 0, 0} = 1 .

x = \frac{4 m ^{2}}{s} .

x = \frac{4 m ^{2}}{s} .

\frac{d}{d x} I_{1} (ϵ; x) I_{2} (ϵ; x) I_{3} (ϵ; x) =

\frac{d}{d x} I_{1} (ϵ; x) I_{2} (ϵ; x) I_{3} (ϵ; x) =

B (x) =

B (x) =

D (x) =

\frac{d}{d x} I_{1 H} (x) I_{2 H} (x) I_{3 H} (x) = B (x) I_{1 H} (x) I_{2 H} (x) I_{3 H} (x),

\frac{d}{d x} I_{1 H} (x) I_{2 H} (x) I_{3 H} (x) = B (x) I_{1 H} (x) I_{2 H} (x) I_{3 H} (x),

G(x)=\left(\begin{array}[]{ccc}H_{1}(x)&J_{1}(x)&I_{1}(x)\\ H_{2}(x)&J_{2}(x)&I_{2}(x)\\ H_{3}(x)&J_{3}(x)&I_{3}(x)\\ \end{array}\right)\qquad\to\qquad\frac{d}{dx}G(x)=B(x)\,G(x)\,.

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

TTP17-021

**Maximal cuts and differential equations for Feynman integrals.

An application to the three-loop massive banana graph

**

Amedeo Primoa,b, 111e-mail: [email protected] and Lorenzo Tancredic, 222e-mail: [email protected]

a* Dipartimento di Fisica ed Astronomia, Università di Padova, Via Marzolo 8, 35131 Padova, Italy *

b* INFN, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy*

c* Institute for Theoretical Particle Physics, KIT, 76128 Karlsruhe, Germany *

We consider the calculation of the master integrals of the three-loop massive banana graph. In the case of equal internal masses, the graph is reduced to three master integrals which satisfy an irreducible system of three coupled linear differential equations. The solution of the system requires finding a $3\times 3$ matrix of homogeneous solutions. We show how the maximal cut can be used to determine all entries of this matrix in terms of products of elliptic integrals of first and second kind of suitable arguments. All independent solutions are found by performing the integration which defines the maximal cut on different contours. Once the homogeneous solution is known, the inhomogeneous solution can be obtained by use of Euler’s variation of constants.

Key words: Feynman graphs, banana graph, sunrise graph, differential equations, unitarity cuts, maximal cut, elliptic integrals

1 Introduction

The study of the mathematical structures that characterize multiloop Feynman integrals has played a crucial role for the most recent developments in the computation of higher order radiative corrections for many processes of direct phenomenological interest at the LHC. The discovery of integration by parts identities [1, 2] together with the differential equations method [3, 4, 5], in particular as extensively applied to compute entire families of massless four-point Feynman integrals [6], allowed to successfully deal with large classes of problems that had remained out of reach before. One step further has been made with the introduction of the concept of a canonical basis [7]. The latter allows to render the computation of multiloop Feynman integrals, which can be expressed in terms of multiple polylogarithms [8, 9, 10], almost entirely algorithmic. In the last years, different methods for the construction of a canonical basis have been proposed, according to the specific properties of the differential equations, such as a linear dependence on the spacetime dimensions [11] or the dependence on a single kinematic variable [12, 13, 14] and concrete steps towards algorithms also valid for multi-scale problems have been made in [15, 16, 17]. The striking simplicity highlighted by the use of canonical bases can be traced back to the rich algebraic structures which can be associated to multiple polylogarithms [18, 19, 20, 21]. More recently, a big effort has been devoted in studying how to possibly extend some of these structures to include more complicated Feynman integrals, which cannot be evaluated in terms of multiple polylogarithms only. These new mathematical functions start to appear at the two-loop order, in particular in the presence of internal massive propagators [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37].

A way to understand the jump in complexity entailed by these integrals is to look at the details of their computation with the differential equations method. From this point of view, Feynman integrals which evaluate to multiple polylogarithms are found to satisfy first order linear differential equations in the external invariants, at least in the limit of even space-time dimensions, $d=2\,n$ , $n\in\mathbb{N}$ 333We recall here that the Laurent coefficients of a Feynman integral expanded around $d=4$ space-time dimensions can be recovered from the corresponding expansion around any even number of dimensions [38]. A way to determine whether the differential equations can be decoupled in $d=2n$ space-time dimensions is described in [39].. What this means is that, it is always possible to write an integral representation (by simple quadrature of the equations) for the Laurent coefficients of their series expansion in $\epsilon=(4-d)/2$ . The situation changes completely when we consider graphs that cannot be expressed in terms of multiple polylogarithms. In this case, the corresponding master integrals are found to satisfy irreducible higher-order differential equations even in the limit $d=4$ , whose most famous example is undoubtedly the second-order differential equation satisfied by the two-loop massive sunrise graph [25].

As it is well known, an integral representation for the solution of a higher-order differential equation can be found only if the full set of its homogeneous solutions is known. Unfortunately, given a generic higher-order differential equation with rational coefficients, except for a limited number of special cases, there is no way to determine the full set of its homogeneous solutions. This constituted for a long time a serious problem in the computation of any Feynman integral of this type. More recently a new picture started to emerge. We showed in [40] that an integral representation for (one of) the homogeneous solutions of the differential equation satisfied by any Feynman integral can be obtained by computing its maximal cut. A similar idea was proposed in [41] in the context of the DRA method. The statement is independent on the number of space-time dimensions $d$ we are interested in, in the sense that the maximal cut computed in any, continuous, number of dimensions $d$ provides an homogeneous solution for its $d$ -dimensional differential equation. Specializing to the physically interesting case of $d=4$ , one obtains at once the homogeneous solution of the corresponding $4$ -dimensional differential equation. It was then shown in [42], that the computation of the maximal cut can be simplified by the use of the so-called Baikov representation [43]. Moreover, the importance of the unitarity cuts for understanding more general mathematical structures of Feynman integrals has been pointed out recently in [44, 45, 46], while in [47] it was shown how to use cuts to derive efficiently the differential equations.

In spite of the increasing effort, until today only a limited number of examples have been considered in the literature. Interestingly, all examples considered could always be reduced to differential equations of at most degree two, whose homogeneous part could always be solved in terms of complete elliptic integrals of first and second kind. Quite in general, for any second-order differential equation, once one homogeneous solution is known, the second can be obtained by a simple quadrature. The same is not true for higher-order equations. Therefore, a very important issue, that was not fully investigated in [40], is whether it is possible to obtain all independent homogeneous solutions from the computation of the maximal cut only. This becomes crucial if one is interested in Feynman integrals which fulfil third- or higher-order differential equations.

In this paper we address this problem and show that the answer to this question is indeed affirmative and that all independent homogeneous solutions can be obtained by evaluating the maximal cut along different independent contours, which do not cross any branch cut of the integrand. We notice here that there is a clear correspondence between this simple idea and the methods described in [48] to count the number of independent master integrals.444Indeed, if a Feynman integral is reduced to $N$ master integrals in $d$ dimensions, the latter will satisfy a system of $N$ coupled differential equations. A complete solution of the latter would therefore require to find $N$ independent homogeneous solutions. While this can be seen very easily already in the case of the two-loop massive sunrise graph, here we move one step forward and consider the first example of a Feynman graph that fulfils an irreducible third-order differential equation: the three-loop massive banana graph. This calculation requires finding three independent homogeneous solutions, for which we show how to derive integral representations from the study of its maximal cut only. We show, moreover, that these integrals can be evaluated explicitly in terms of products of elliptic integrals of first and second kind of complicated arguments. This result, which is very non-trivial, was expected. In fact, the third-order differential equation satisfied by this graph has been studied long ago by G.S. Joyce in the context of cubic lattice Green functions [49], where it was shown that this equation is a so-called symmetric square and its solution can be therefore expressed as products of solutions of a lower second-order differential equation. Finally, we recall the reader that an alternative calculation of this graph has been presented in [50], where it was shown that the three-loop banana graph in $d=2$ can be written as an elliptic three-logarithm.

The rest of the paper is organized as follows. In Section 2 we reconsider in detail the well known case of the two-loop massive sunrise graph, showing how the two independent homogeneous solutions of its second-order differential equation can be found by integrating the maximal cut over the only two independent contours that one can build without crossing any branch cut of the integrand. We switch then to the more interesting case of the three-loop massive banana graph in Section 3, where we derive the system of three-coupled differential equations that it satisfies. In Section 4 we focus on the system of differential equations around $d=2$ space-time dimensions, which is known to be equivalent to the corresponding expansion close to $d=4$ [38]. We compute the maximal cut of the three-loop banana graph in Section 5 and show that by evaluating it along different contours we obtain at once all independent solutions of the homogeneous differential equations. In Section 6 we show that the same third-order differential equation is a symmetric square, we describe its solution as derived in [49] and discuss its equivalence to ours. Finally in Section 7 we put everything together and compute the inhomogeneous solution, continuing it analytically to the whole phase space. We also provide some appendices where we describe some of the technical calculations required in the main text.

Note Added

Right before the completion of this work an interesting paper appeared on the arXiv [51], containing similar findings to the one presented here. In [51] it is shown, with many explicit examples, how to build the full set of independent homogeneous solutions of the differential equations by integrating the maximal cut on different regions, for generic values of the space-time dimensions $d$ . This is done efficiently using Baikov representation. While our conclusions are similar, in our paper we do not dwell on general methods for an efficient computation of the cut in any numbers of dimensions, as we are primarily interested in the calculation of the Laurent coefficients of the master integrals close to $d=4$ and the complexity of this calculation is determined entirely by the value of the homogeneous solutions at $d=2n$ , $n\in\mathbb{N}$ . Instead, we apply for the first time these ideas to solve a system of three coupled differential equations, which require a generalization of the methods used for the two-loop massive sunrise graph and other similar examples. From this point of view, we believe our results are complementary to the ones presented in [51].

2 Revisiting the two-loop massive sunrise graph

Before looking in detail at the case of the three-loop massive banana graph, it is useful to reconsider the better understood case of the two-loop massive sunrise. This will help us to illustrate some of the techniques used later for the three-loop banana graph in a simpler environment.

As it is well known, the sunrise graph satisfies a coupled system of two linear differential equations [25]; its homogeneous solution is therefore a $2\times 2$ matrix, whose entries are linear combinations of square-roots and complete elliptic integrals of first and second kind, see for example [34]. The entries of this matrix were determined for the first time in [25] by studying the imaginary part of the graph. Taking inspiration from this calculation, in [40] we showed that the computation of the maximal cut of any family of Feynman integrals naturally generalizes the approach of [25] and allows one to determine an integral representation for an homogeneous solution of its system of differential equations irrespective of its degree. As argued in [40], if an homogeneous solution can be determined in terms of complete elliptic integrals, it is then straightforward to find the second solution using the properties of the elliptic integrals. While this is very useful, it cannot be straightforwardly extended to more general cases which involve higher-order differential equations and are not expected to be solved in terms of elliptic integrals only.

The goal of this section is to show how the study of the maximal cut of the sunrise graph provides at once both independent solutions. While the validity of this claim is somehow obvious, its consequences are non-trivial and far-reaching. This remains valid, in fact, for any Feynman integral and provides therefore an effective method for constructing integral representations for all independent homogeneous solutions of the differential equations satisfied by the graph, regardless of their complexity. In particular, we will show explicitly how this is done by integrating the maximal cut of the sunrise graph along the only two linearly independent integration contours which do not cross any singularity of the integrand.

Let us consider the equal-mass sunrise graph

[TABLE]

where $u=p^{2}/m^{2}$ . The graph in $d=2$ space-time dimensions satisfies the following second-order differential equation

[TABLE]

where we neglected the inhomogeneous terms which are irrelevant here and we set $S(2;u)=S(u)$ . As it is well known, the maximal cut of the sunrise graph in $d=2$ can be written as

[TABLE]

where we use the notation ${\rm Cut}(S(u))$ for the maximal cut of $S(u)$ and we have not fully specified neither the integration contour $\mathcal{C}$ nor the sign of the argument of the root. We claim that the integration along any contour $\mathcal{C}$ which does not cross any branching point of the integrand produces a solution of (2.2). In particular, we will see that there are only two possible independent contours of such type and that by integrating along them we get at once both independent solutions of (2.2).

First of all, the square-root has four branching points. By choosing $u>9$ we have

[TABLE]

The ordering of the branching points depends on the value of $u$ , but the argument used below does not depend on it. Given the four branching points it should be obvious that, depending on the sign that we pick in (2.3), there are four possible integration contours which we can draw without crossing the branch cuts. If we choose the plus sign, the integrand develops a branch cut for $0<b<4$ and $(\sqrt{u}-1)^{2}<b<(\sqrt{u}+1)^{2}$ . If we pick the minus sign the branches are for $-\infty<b<0$ , $4<b<(\sqrt{u}-1)^{2}$ and $(\sqrt{u}+1)^{2}<b<+\infty$ . In the first case, i.e. picking a plus sign, we can clearly draw the two contours $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ depicted in Figure 1a. The third contour, $\mathcal{C}_{\infty}$ , is instead equivalent to the sum of the two, and we will need it later on. In the second case, we can draw instead only one single contour, see Figure 1b, giving a total of three apparently different possibilities.

Nevertheless, it is easy to see [25] that two of the three contours are equivalent. Let us consider the following integral

[TABLE]

where $\mathcal{C}_{\infty}$ is the contour at infinity encircling the four roots (2.4) depicted in Figure 1a. Clearly the integral is zero since the contour does not contain any poles and the integrand goes as $1/b^{2}$ when $b\to\infty$ . On the other hand, with the sign of the root in (2.5), the integrand has two cuts, one for $0<b<4$ and the other for $(\sqrt{u}-1)^{2}<b<(\sqrt{u}+1)^{2}$ . We can then imagine to shrink $\mathcal{C}_{\infty}$ to encircle the two branch cuts and we get from (2.5)

[TABLE]

which proves that the two integrals are not independent. In this way we are left with two independent contour integrals which provide precisely the two independent solutions of (2.2), say $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ . Now, by shrinking $\mathcal{C}_{1}$ , $\mathcal{C}_{2}$ and $\mathcal{C}_{3}$ on the corresponding branch cuts one finds an equivalent representation as one-dimensional real integrals

[TABLE]

where the sign of the roots on the right hand side of the formulas is chosen to deal with real integrals. Alternatively, the contour $\mathcal{C}_{2}$ can also be sent to infinity providing a second integral representation for the integral

[TABLE]

We have therefore determined two independent solutions as the integrals over the two independent contours, say $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ , which in turn we can written as one-dimensional real integrals as follows

[TABLE]

It was shown in [25, 34] that these functions are indeed the two independent solutions of (2.2).

Clearly, the analysis we performed is completely general and does not depend on the precise position of the branching points (2.4). In fact, given a generic integral of the form

[TABLE]

where $a_{1}$ , …, $a_{4}$ are four distinct roots such that

[TABLE]

the same considerations apply and one finds that there are only two independent contours which are equivalent to the following one dimensional real integrals

[TABLE]

As for the explicit case of the sunrise with equal masses, we chose the sign in the root in order to deal with real integrals everywhere.

As a last remark, the integrals in Eqs. (2.11) are nothing but complete elliptic integrals of the first kind. To see this, we perform the two standard changes of variables for the two integrals respectively

[TABLE]

and obtain

[TABLE]

where

[TABLE]

is the elliptic integral of the first kind and

[TABLE]

Indeed, a standard result of the theory of the complete elliptic integrals shows that $K(w)$ and $K(1-w)$ satisfy the same second-order differential equation, of which they constitute the two independent solutions.

The analysis carried out in this section might seem somewhat redundant, as the theory of the elliptic integrals has been very well understood for a long time. Nevertheless, when considering the three-loop banana graph, we will see that many of the ideas and of the results derived here can be directly borrowed or trivially extended to more complicated cases. We believe that this will make our analysis in this much less trivial case, much more transparent.

3 The three-loop massive banana graph

We consider the three-loop two-point integral family defined by

[TABLE]

with $p^{2}=s\neq 0$ . Our integration measure is defined as

[TABLE]

such that the one-loop tadpole integral reads

[TABLE]

The family of integrals in (3.1) with $a_{5},...,a_{9}<0$ is often referred to as three-loop banana graph or three-loop sunrise graph. A first step towards the computation of these integrals is to use integration by parts identities to reduce them to a basis of master integrals. In this case, this step can be very easily achieved with any of the publicly available reduction codes like, for example, Reduze 2 [52]. There is only one sub-topology, the three-loop massive tadpole, which is reduced to one single master integral. The top-topology, instead, in the limit of equal internal masses can be reduced to three independent master integrals. The choice of the master integrals is of course arbitrary. Since massive two-point functions are IR finite, they are best studied in the vicinity of $d=2$ instead of $d=4$ , since for $d=2$ the scalar amplitude is also UV finite. We define therefore $\epsilon=(2-d)/2$ and we pick the following basis

[TABLE]

where the normalization factors have been chosen for later convenience. For the three-loop tadpole we chose instead the master integral $I_{2,2,2,0,0,0,0,0,0}$ , which in our normalization becomes simply

[TABLE]

The dependence of the master integrals on the momentum transfer $p^{2}=s$ can be conveniently parametrized in terms of the dimensionless variable

[TABLE]

The system of differential equations in $x$ satisfied by $\mathcal{I}_{1}$ , $\mathcal{I}_{2}$ and $\mathcal{I}_{3}$ is given by

[TABLE]

where $B(x)$ and $D(x)$ are $3\times 3$ matrices that do not depend on $\epsilon$ ,

[TABLE]

and the inhomogeneous term is proportional to the massive tadpole (3.5). We stress once more that we study the solution of the differential equations in the vicinity of $d=2$ and in our conventions $\epsilon=(2-d)/2\,.$

As it is easy to see, the structure of the system of differential equations is strikingly simple, as it is characterized by only four regular singular points $x=0$ , $x=1/4$ , $x=1$ and $x=\pm\infty$ , which correspond, respectively, to $s=\pm\infty$ , $s=16m^{2}$ , $s=4m^{2}$ and $s=0$ . The structure of the singularities resembles closely that of the simpler two-loop massive sunrise graph, see for example [25]. Given the simplicity of the equations, it is particularly interesting to investigate the class of functions that are needed for their solution and in which sense these functions generalize the ones required for the integration of the two-loop sunrise graph. In the two-loop case, the homogeneous system of differential equations admits as solutions complete elliptic integrals of the first and second kind. Indeed, as we will see in the next section, the three-loop case can be solved in terms of products of complete elliptic integrals of first and second kind.

4 The homogeneous system

Before embarking in the solution of the differential equations, let us first recall which ingredients are needed to solve a $3\times 3$ coupled system, as the one in Eq. (3.7). Since we are interested in calculating the solution as a Laurent expansion in $\epsilon$ , the first step consists in solving the homogeneous system for $\epsilon=0$

[TABLE]

where the suffix $H$ indicates that we are limiting ourselves to its homogeneous piece. The solution of Eq. (4.1) requires finding a $3\times 3$ matrix, say $G(x)$ , defined such that

[TABLE]

The inverse of the matrix $G(x)$ is

[TABLE]

where $W(x)=\det\left(G(x)\right)$ is the Wronksian of the system (4.1). As it is well known $W(x)$ satisfies the first order differential equation

[TABLE]

which is often referred to as Abel’s identity. From the matrix $B(x)$ defined in Eq. (3.11) we get at once

[TABLE]

which admits the very simple solution

[TABLE]

where we assumed for simplicity to work in the Euclidean region, $x<0$ . The value of the integration constant $c_{1}$ can be determined only once the exact form of the solutions is known.

The simplicity of the Wronskian (4.9) is a general feature and can be easily understood, irrespective of the complexity of the matrix $G(x)$ , from the fact that it always satisfies a first order differential equation (4.7). This has interesting consequences. Quite in general, let us consider the homogeneous system Eq. (4.1) and write it in components as

[TABLE]

Note that we have written the system for a generic set of functions $f_{i}(x)$ , since all considerations made here can be applied to any other $n\times n$ system of differential equations and are not limited to the particular set of master integrals $\mathcal{I}_{i}(x)$ chosen above. Let us limit ourselves to the diagonal part of the system, which reads

[TABLE]

It is straightforward to solve these equations by quadrature as

[TABLE]

where $x_{0}$ is an irrelevant boundary condition. With this we can define a new basis of master integrals

[TABLE]

which fulfil a new system of differential equations

[TABLE]

whose matrix $\widetilde{B}(x)$ is now traceless ${\rm Tr}\left(\widetilde{B}(x)\right)=0$ by construction. Due to Abel’s identity (4.7), the Wronskian $W(x)$ of the new system becomes now particularly simple

[TABLE]

This implies in turn that, irrespective of the complexity of the functions appearing in the solution matrix $G(x)$ (4.2), there always exists (at least) one trivial combination of them such that

[TABLE]

In the well known case of elliptic integrals, which satisfy second-order differential equations, this relation reduces to the so-called Legendre relation.

The determination of the entries of the matrix (4.2) is in general very non-trivial. We will show how to do this using the information coming from the maximal cut in the next section, as it was first suggested in [40]. Assuming that the matrix $G(x)$ is known, one can define a new basis of master integrals $M_{i}$

[TABLE]

which by construction fulfil the system

[TABLE]

Once the system is in form (4.18), its solution order by order in $\epsilon$ reduces, at least in principle, to a quadrature. In particular, the matrix $G^{-1}(x)D(x)G(x)$ (combined with the integral over the inhomogeneous term in (4.18)) contains all information concerning the class of functions needed to iterate the solution at every order in $\epsilon$ . The problem becomes then that of classifying these functions and of understanding their analytic and algebraic properties. Indeed, for a general problem it will not be possible to solve these integrals in terms of known special functions. In our case, nevertheless, we managed to write the matrix $G(x)$ for the three-loop banana graph in terms of relatively simple products of elliptic integrals and it would be interesting to investigate whether a generalization of the elliptic polylogarithms introduced to represent the two-loop massive sunrise could be used also in this case, as the calculation in [50] seems to imply.

5 The maximal cut and the homogeneous solution

Our goal is now to solve the homogeneous system (4.1) and determine analytically the form of the matrix $G(x)$ . A $3\times 3$ coupled system can always be rephrased as a third-order differential equation for any of the master integrals. It is natural to do this for the scalar amplitude $\mathcal{I}_{1}(\epsilon;s)$ , whose third-order homogeneous differential equation in $d=2$ space-time dimensions reads

[TABLE]

The solution of this equation in terms of three linearly independent functions provides the first row of the matrix $G(x)$ , see Eq. (4.2), while the remaining two rows can by obtained by differentiation with respect to $x$

[TABLE]

Solving Eq. (5.1) is non-trivial. This differential equation has been studied long ago in the context of the calculation of cubic lattice Green functions [49], where it was solved in terms of products of two elliptic integrals. In order to re-derive this result we will use the information coming from the maximal cut. As it was shown in [40], given a set of master integrals, their maximal cut singles out by construction the homogeneous part of the differential equations. We use the notation ${\rm Cut}(\mathcal{I}_{j}(x))$ for the maximal cut of $\mathcal{I}_{j}(\epsilon;s)$ evaluated in $\epsilon=0$ and obtain at once

[TABLE]

or equivalently

[TABLE]

Since Eq. (5.5) admits three independent solutions, and all of them are required in order to construct the solution of the system, the really interesting question becomes how one can obtain all of them from the maximal cut of $\mathcal{I}_{1}(x)$ only. As we showed in Section 2 for the simpler case of the two-loop sunrise graph, this can be done by integrating the maximal cut along the independent contours that do not cross any branch cuts of the integrand. Let us see how this works for the present case.

First of all, given the definition of the banana graph in Eq. (3.1), we note that the scalar amplitude in $d=2$ dimensions is finite and can be written as an integral over two one-loop bubbles

[TABLE]

In this way, the computation of the maximal cut for the three-loop banana graph can be greatly simplified. We start with the massive one-loop bubble in $d=2$ dimensions

[TABLE]

Its maximal cut can be easily computed to be, up to an irrelevant overall normalization constant,

[TABLE]

Using Eq. (5.8) the maximal cut for the banana graph can then be written as

[TABLE]

where we will specify the different choices for the contour $\mathcal{C}$ later on. In $d=2$ a way to simplify this integral is to parametrize the loop momentum and the external one in terms of two massless momenta $p_{1}^{\mu}$ and $p_{2}^{\mu}$ as follows

[TABLE]

With this parametrization we get

[TABLE]

and where $\mathcal{C}^{\prime}$ represents now a still unspecified two-dimensional contour in the four-dimensional complex hyperplane spanned by the variables $a,b$ . In this way the integral becomes, up to a multiplicative constant

[TABLE]

where we used $x=4\,m^{2}/s$ and introduced the polynomial

[TABLE]

Determining explicitly which integration contours provide the independent solutions is a mathematically interesting problem, related to the dimension of the so-called cohomology group associated to the variety described by the curve $R(a,b,c)=0$ . Instead of embarking into complicated mathematical considerations that go beyond the scope of this paper, we can try and see what we can say using simple mathematics. Let us consider integral (5.12). If we allow $a$ and $b$ to assume complex values, we are effectively integrating on a complex two-dimensional hypersurface embedded in the four-dimensional complex space spanned by $a$ and $b$ . There are two fundamental questions we should answer. The first is, indeed, which contours provide a solution of the differential equation (5.1). Once we have determined all possible contours, the second question is which ones are linearly independent and provide therefore independent solutions. The answer to the first question is quite general, as we claimed in Section 2. Any complex contour that does not cross any branch cuts will do the job. The second question, instead, has to do with recognizing which of the allowed contours are linearly independent from each other.

It is of course very difficult to picture a two-dimensional surface embedded in a four-dimensional space. Instead, it is useful to focus on the real two-dimensional plane spanned by the real coordinates associated to $(a,b)$ . It is trivial to see that the square-root in Eq. (5.12) has zeros for $a=0$ , $b=0$ , $a=-1$ , $b=-1$ , $a=x/b$ and $a=x/(b+1)-1$ . We draw these sets of points in Figure 2 as continuous lines of different colors. Similarly to the simpler case of the two-loop sunrise studied in Section 2, the square root in Eq. (5.12) can develop a branch cut every time that one crosses one of these lines. Now it is easy to convince oneself that a two-dimensional complex contour that does not cross any branch cuts of the integrand cannot cross in particular any of the lines in Figure 2. We can then imagine to shrink these contours to different two-dimensional regions in the real plane bounded by the branches drawn in Figure 2. Integrating the maximal cut (5.12) in any of these regions will provide us with a viable solution of (5.1).

Now that we know which contours are allowed, we should determine which ones of them are independent and provide us therefore with the independent solutions of (5.1). This analysis is simplified by considering integral (5.12) as an iterated integral over two one-dimensional contours in $da$ and $db$ 555Notice that the starting integral is symmetric in $a$ and $b$ , so the analysis that follows can be equally well repeated inverting the role of the two variables.

[TABLE]

where we introduced the abbreviations

[TABLE]

From Eq. (5.14) we see that the integral in $da$ is an integral over a square root of a quartic polynomial, which defines an elliptic curve. We studied this integral in section 2. That analysis implies that for every fixed value of $x$ and $b$ (which in turn determine an ordering of the roots (5.15)) there are two independent contours $\mathcal{C}_{a}$ which correspond to the two periods of the elliptic curve, see Eqs (2.11). Without any loss of generality, we assume for the time being $1/2<x<1$ , i.e. $4m^{2}<s<8m^{2}$ . This choice is required for the analysis below, but the result obtained is of course independent of it and can be analytically continued to any other value of $s$ .

Given this choice, we start by dividing the $(a,b)$ -plane into five different regions depending on the value of the variable $b$

[TABLE]

This is useful since in each region the branching points $a_{j}$ given in Eq. (5.15) are ordered differently and it allows us to define real building blocks to construct our solutions. This is done as follows. First of all, for each of the regions (5.16), we can define two independent functions from the two corresponding contours $\mathcal{C}_{a}$ identified in (2.11), and we get a total of $10$ possibilities

[TABLE]

where the sign in the square-root is chosen to deal always with real functions. It is very important to realize that these integral representations are not unique. Indeed, as we discussed at length in Section 2, there are different equivalent representations for the integrals in $da$ , see (2.11). This will turn out to be crucial for the considerations below.

We should notice that many of these functions are related to each other by simple symmetry transformations. First of all, the integrand in Eq. (5.12) is symmetric under $\{a\leftrightarrow b\}$ and $\{a\to-a-1,\;b\to-b-1\}$ (and of course under a combination of the two symmetries). It is immediate to see that under $\{a\to-a-1,\;b\to-b-1\}$ region I and V and region II and IV are related to each other and we find

[TABLE]

The effect of the symmetry under $\{a\leftrightarrow b\}$ is less immediate. Let us take for example $f_{1}^{\rm V}(x)$ . Using $R(a,b,x)=R(b,a,x)$ , we rename $a$ and $b$ and exchange the order of integration obtaining

[TABLE]

where in the last step we used (5.22). With a similar calculation one can show that

[TABLE]

leaving in this way four independent functions, which we can choose to be

[TABLE]

The four functions (5.25) are not all solutions of the third-order differential equation. We should remember, in fact, that a solution is obtained when we integrate Eq. (5.12) in a region of the $(a,b)$ -plane bounded by the branch cuts. It is easy to see from Figure 2 and from the definition of the functions (5.20, 5.21) that while $f_{1}^{\rm V}(x)$ , $f_{2}^{\rm V}(x)$ indeed fulfil this requirement, $f_{1}^{\rm IV}(x)$ and $f_{2}^{\rm IV}(x)$ apparently do not. This could constitute a problem since, in order to solve a third-order differential equation we need three independent solutions, while this argument seems to suggest we can find only two. Indeed, there is a subtlety. To understand it, let us look more closely at Figure 2. There, the integrals corresponding to these four functions are depicted as shaded areas of different colors. One should remember that, for any given region in $b$ , different choices of integration boundaries in $a$ produce equivalent results, as summarized in (2.11). In particular, focussing on regions $\rm IV$ and $\rm V$ , we see immediately that in both cases the integrals corresponding to the areas $1$ and $3$ are actually identical. In the language of the functions defined above, this means for example that

[TABLE]

and the identity follows from Eq. (2.11), which specialized in this case becomes simply

[TABLE]

As it is clear from the figure, the second integral representation in (5.26) is now integrated in a region bounded by the branch cuts, which means that $f_{1}^{\rm IV}(x)$ actually also fulfills our requirements and it is therefore a solution of the third-order differential equations. The same cannot be said $f_{2}^{\rm IV}(x)$ , which cannot be rewritten as an integral over a region of plane bounded by branch cuts only. The complete set of three independent solutions is therefore given by $f_{1}^{\rm V}(x)$ , $f_{2}^{\rm V}(x)$ and $f_{1}^{\rm IV}(x)$ .

Interestingly, if we want to build a viable solution of the equation using the function $f_{2}^{\rm IV}(x)$ , we should consider the combination

[TABLE]

As this function is defined by integrating in a region delimited by the branch cuts of the root, it must be a solution of the third-order differential equation. On the other hand, since we have already determined the three independent solutions, it must be possible to write it as a linear combination of $f_{1}^{\rm V}(x)$ , $f_{2}^{\rm V}(x)$ and $f_{1}^{\rm IV}(x)$ . Indeed, using (5.22), (5.23) and (5.24) one finds

[TABLE]

showing that the solution is indeed not linearly independent. Of course until now we have not proved directly that the functions determined above actually solve the third order differential equation. We will do it later on once we have found a more convenient representation for them.

5.1 A basis of unit leading singularity

Before embarking in the explicit evaluation of the integrals above, let us pause to interpret our results. In the previous section we have showed that the maximal cut of the three-loop banana graph can be evaluated on three independent contours which do not cross branch cuts of the integrand and we claimed that the latter provide the three independent solutions of its differential equation. This provided us with all ingredients to build the complete homogeneous solutions of the original $3\times 3$ system (4.1), as a matrix $G(x)$ (4.2). As we showed in Eq. (4.17), the matrix $G(x)$ allows us to rotate the original basis of master integrals $\mathcal{I}_{1}(x)$ , $\mathcal{I}_{2}(x)$ and $\mathcal{I}_{3}(x)$ , onto a new basis $M_{1}(x)$ , $M_{2}(x)$ and $M_{3}(x)$ which satisfies differential equations almost entirely $\epsilon$ -factorized.

The new basis defined in this way, presents a remarkable (but obvious) property. It’s maximal cut, computed along the three independent integration contours, is unity (or zero) by definition. This can be easily seen as follows. Inverting Eq. (4.17) we have

[TABLE]

Now, from the results of the previous section, we can identify

[TABLE]

where we used the notation ${\rm Cut}_{\mathcal{C}}(\mathcal{I}_{j}(x))$ for the maximal cut of $\mathcal{I}_{j}(x)$ computed along the contour $\mathcal{C}$ , which gives rise to the integration over the corresponding regions discussed above. The remaining functions can be obtained by differentiating the ones above as in (5.2, 5.3), or alternatively by computing the maximal cuts of the other two master integrals along the same three contours

[TABLE]

Using now Eq. (5.29, 5.30, 5.31) we find at once

[TABLE]

This result is important. We can imagine, in fact, to associate to any family of master integrals which fulfil a set of $n$ (in our case 3) coupled differential equations, a $n\times n$ matrix whose entries are given by the maximal cut of the integrals evaluated along the $n$ independent integration contours. For a generic basis of master integrals, this is by definition the matrix $G(x)$ which solves the homogeneous system of differential equations. For the rotated basis $M_{j}(x)$ , Eq. (5.32) shows that remarkably this “matrix of maximal cuts” reduces to the identity matrix.

With this, we can therefore imagine a way to generalize the idea originally presented in [7], where it was claimed that master integrals with unit leading singularity are expected to fulfil canonical differential equations. Suppose we are considering a basis of $n$ master integrals which fulfil $n$ coupled differential equations. In this case, a basis which fulfils $\epsilon$ -factorized differential equations must have unit leading singularity in the sense above, i.e. the matrix which contains as entries the maximal cut of the master integrals evaluated along all independent integration contours must be equal to the identity matrix. Indeed, it becomes clear that, from a practical point of view, there is no real difference between finding a basis of unit leading singularity and actually solving the homogeneous system of differential equations.

5.2 The homogeneous solutions as product of elliptic integrals

We go back now to the explicit form of the homogeneous solution for the three-loop banana graph. The analysis above has allowed us to determine the three independent solutions in form of two-fold integral representations. We might now ask ourselves if these integrals can be performed in terms of known functions. The answer is indeed affirmative.

First of all, in order to proceed, it is useful to perform explicitly the integration in $da$ and obtain a one-fold integral representation for the solutions. Using the change of variables given in (2.12) on the four function (5.25) and after a bit of algebra we find respectively

[TABLE]

The integral representations (5.33, 5.34, 5.35) are already more convenient for numerical integration and analytical manipulations. Nevertheless we can do better and compute the three integrals explicitly in terms of products of elliptic integrals of the first kind only. The manipulations are non-trivial and we found easier to show how this is done on the last two integrals, $f_{1}^{\rm V}(x)$ and $f_{2}^{\rm V}(x)$ . Similar manipulations should be possible also for the first functions, $f_{1}^{\rm IV}(x)$ , but as it will become clear, we will not need to perform them directly. Let us then consider Eqs. (5.33, 5.34) and perform the change of variables

[TABLE]

such that the two integrals become

[TABLE]

Eqs. (5.37,5.38) can be put in standard form by introducing three parameters $\alpha$ , $\beta$ and $\gamma$ defined such that

[TABLE]

In general, for a given value of $x$ , four different pairs of solutions exist to these equations, and we can choose any of those for what follows. For definiteness, we pick

[TABLE]

where we are assuming $0<x<1$ . For any given pair of solutions, the integrals read

[TABLE]

Integrals (5.41, 5.42) are now in standard form. In particular, the calculation of (5.42) is discussed in [53], see Eq.(33) therein. Suitable extensions of the methods described there allow us to compute also also integral (5.41). The results read

[TABLE]

where we have defined

[TABLE]

It is easy to prove by direct calculation that Eqs. (5.43, 5.44) solve the third order differential equation (5.1), for every choice $\alpha,\beta$ and, in particular, for the choice we made in (5.40). Even if the results for the integrals Eqs. (5.41, 5.42) can be found in standard tables of integrals, their calculation remains in our opinion not entirely straightforward and we report it therefore in Appendix B following closely the methods described in [53].

Even if we have computed only two of the solutions, by direct inspection of Eqs. (5.43, 5.44), we can easily identify the three independent solutions of the third-order differential equation (5.1) as the three functions

[TABLE]

with $k_{\pm}$ given by (5.45) together with (5.39) and $x=4m^{2}/s$ . These results allow us to fix completely the first row of the matrix of solutions (4.2). The other two rows can then be obtained by Eqs. (5.2), (5.3). We do not report the results here for brevity, as we will use a different and more compact representation later on. We can nevertheless use these solutions to compute the Wronskian and find, as expected from Abel’s formula (4.7)

[TABLE]

where of course the overall normalization constant depends on the explicit normalization choice made on Eqs. (5.46). Inspecting our three solutions, it is natural to wonder what would happen considering one further combination

[TABLE]

It is simple to prove by direct calculation that also Eq. (5.48) solves the differential equation (5.1). Of course, since a third-order differential equation admits only three independent solutions, this last solution cannot be linearly independent from the previous three. Indeed, it is easy to check numerically (for example using the PSLQ algorithm) that this is true. Since the four functions develop imaginary parts for $x<0$ or $x>1/4$ , the exact relation between the solutions depends on the value of the variable $x$ and on the convention picked for their analytic continuation. With the choice (5.40) and taking $0<x<\frac{1}{4}$ ( $s>16m^{2}$ ), for which all solutions Eqs. (5.46,5.48) are real valued, one finds

[TABLE]

One last comment is in order. We have computed explicitly two out of the three functions that we claimed constitute the three solutions of the third-order differential equation. We have been lucky enough to be able to extract from these two functions a representation for all three independent solutions of the equation as products of elliptic integrals (5.46). This of course has to be considered as an accident. One might wonder what would the calculation of the remaining function $f_{1}^{\rm IV}(x)$ have produced. Indeed, if it is true that it is also a solution of the equation, and if the solutions chosen in (5.46) are really independent, it must be possible to write $f_{1}^{\rm IV}(x)$ as linear combination of the latter. Instead of performing the integration explicitly, it is simple to prove this by use of the PSLQ algorithm. As before, the relation depends on the value of $x$ . The function $f_{1}^{\rm IV}(x)$ is real-valued for $1/4<x<1$ , while the three solutions (5.46) are real-valued for $0<x<1/4$ . Picking then $1/4<x<1$ and assuming $x\to x+i\,0^{+}$ , one finds easily the following relation

[TABLE]

showing that, as expected, also $f_{1}^{\rm IV}(x)$ is a solution of the third-order differential equation satisfied by the banana graph.

6 The third-order differential equation as a symmetric square

In the previous section we have showed that different solutions of the third-order differential equation (5.1) can be found by integrating the maximal cut of the three-loop massive banana graph along independent integration contours. Here we want to elucidate the relation of the solutions found above with an alternative set of solutions found by G.S. Joyce [49]. For simplicity, let us rewrite here the third-order homogeneous differential equation satisfied by the three-loop banana graph

[TABLE]

This equation has a remarkable property, i.e. it is a so-called symmetric square. Completely in general, let $L_{3}(x)$ be a third-order differential operator

[TABLE]

The operator $L_{3}(x)$ is a symmetric square if its three independent solutions can be written as

[TABLE]

where the two functions $f_{1}(x)$ and $f_{2}(x)$ are in turn solutions of a second-order differential operator

[TABLE]

We have therefore

[TABLE]

and correspondingly

[TABLE]

Testing whether a third-order differential operator is a symmetric square and building the corresponding second-order differential operator is very simple. Starting from the coefficients of $L_{3}(x)$ in (6.2) we can build the following two combinations

[TABLE]

where $\alpha_{1}^{\prime}(x)=d\alpha_{1}(x)/dx$ . Now, if the following relation is satisfied

[TABLE]

then the differential operator $L_{3}(x)$ in (6.2) is the symmetric square of the corresponding second-order differential operator

[TABLE]

It is straightforward to check whether our third-order differential equation (5.1) is a symmetric square. Starting with

[TABLE]

we immediately find the corresponding coefficients

[TABLE]

and indeed one can verify that

[TABLE]

This shows that the solutions of the third-order differential equation satisfied by the three-loop banana graph (5.1) can be written as symmetric square combinations of the two solutions of the following second-order differential equation

[TABLE]

which can be then solved in terms of a class of special functions called Heun functions [49]. Interestingly enough, one can show that such Heun functions can be rewritten as a product of elliptic integrals of suitable arguments. The result is very non-trivial and we refer to [49] and references therein for details. The three solutions found there read

[TABLE]

where we defined

[TABLE]

Moreover, in [49] it is shown that a relation exists between the elliptic integrals above. If we assume $x>1$ than the functions (6.14) are explicitly real and the relation reads

[TABLE]

This representation of the solutions is somewhat more compact than what we found above in (5.46), but of course the two sets of solutions must be equivalent since they solve the same third-order differential equation. Proving by algebraic manipulations that the two set of solutions can be written one in terms of the other is non trivial, in particular since the exact relations among the two depends on the region in $x$ which we consider. Indeed, every time one crosses a branching point of the differential equation (5.1), each of the three solutions gets in general mapped to a linear combination of the same three functions. To give an example, we consider again the region $x>1$ where the functions (6.14) are real. The solutions found in (5.46) instead are complex and we need a prescription, which we choose for definiteness to be $x\to x+i0^{+}$ . By using PSLQ one then finds (with virtually arbitrary precision) that the following relations are satisfied

[TABLE]

showing that the two sets of solutions are indeed equivalent. The remaining two columns of the matrix of solutions $G(x)$ can then be determined differentiating (5.46) as in Eqs. (5.2, 5.3) . Using these solutions one find for the Wronskian

[TABLE]

We obtained in this way two equivalent representations for the entries of the matrix of homogeneous solutions $G(x)$ defined in (4.2). Still we are not quite done. Since the differential equation (5.1) has four singular points $x=0$ , $x=1/4$ , $x=1$ and $x=\infty$ , its solutions can develop branch cuts crossing any of these points. It is easy to see that the solutions built in this section are real for $x>1$ but they develop an imaginary part whenever $x<1$ . Moreover they have discontinuities in all other singular points. On the other hand, the solutions found in (5.46) turn out to be real only for $0<x<1/4$ . In order to properly analytically continue the results for every value of $x$ we will need to build other solutions, similar to (6.14) or (5.46) , but which are real in the remaining regions $(-\infty,0)$ , $(1/4,1)$ . Indeed, in every region the solutions can be built taking simple linear combinations of (5.46) or (6.14). Many different combinations are possible and we refer to Appendix A for details on one possible choice which we found convenient. To indicate these different sets of solutions we introduce the notation $H_{k}^{(a,b)}(x)$ , $J_{k}^{(a,b)}(x)$ , $I_{k}^{(a,b)}(x)$ , where the superscript indicates that the corresponding solution is real for $a<x<b$ . The corresponding matrix of solutions will then be indicated as $G^{(a,b)}(x)$ . We normalize our solutions for all values of $x$ in such a way that the corresponding Wronskian $W^{(a,b)}(x)$ is equal to (6.18), up to a possible overall factor $i$ which is required when the argument of the square-root becomes negative.

7 The inhomogeneous solution

In this section we make use of the homogeneous solutions of the system of differential equations (3.7), which have been studied in Sections 5-6 and collected in Appendix A, in order to build the inhomogeneous solution as a series expansion around $d=2$ .

The following discussion holds for any of the kinematic regions $a<x<b$ located by the four singular points of the differential equations and we will keep giving as understood the superscript $(a,b)$ .

Both basis of master integrals $\mathcal{I}_{i}(x)$ and $M_{i}(x)$ are finite in $d=2$ and they can be Taylor-expanded as

[TABLE]

By substituting Eq. (7.1) into (4.18), we obtain a particularly simple set of chained first order differential equations for the coefficients $M_{i}^{(k)}(x)$ , which reads

[TABLE]

and

[TABLE]

From Eq.(7.3) we see that the solution for the master integrals $M_{i}(x)$ has a manifest iterative structure, since each coefficient $M_{i}^{(n)}(x)$ can be simply written as an integral of the lower order terms $M_{i}^{(k-1)}(x)$ , convoluted with the integration kernel $G^{-1}(x)D(x)G(x)$ . Together with the integration of the lowest order (7.2), this kernel specifies the class of functions required at every order in $\epsilon$ .

Once $M_{i}^{(n)}(x)$ have been determined, the corresponding term of the $\epsilon$ -expansion of the original master integrals $\mathcal{I}_{i}(x)$ can be obtained by applying the rotation matrix $G(x)$ back to the integrals $M_{i}(x)$ , according to the definition (4.17),

[TABLE]

In the remaining of this section, we will limit ourselves to the determination of the order zero terms $\mathcal{I}_{i}^{(0)}(x)$ . The latter, by definition, are not the entire story, as their calculation does not require to integrate over the kernel $G^{-1}(x)D(x)G(x)$ . While we could extend the methods described in this section in order to provide integral representations also for the higher orders, we do not find this particularly useful and postpone this problem to later work. A complete solution of this problem, in fact, would require to understand and classify the properties of the functions defined by repeated integrations over the kernel above.

The differential equations (7.2) can be readily solve solved by quadrature, producing

[TABLE]

where the integration base-point $x_{0}$ can be arbitrarily chosen and the integration constants $c_{i}^{(0)}$ have to be fixed by imposing suitable boundary conditions. The integrands $\mathcal{R}_{i}(x)$ are combinations of products of two homogeneous solutions which originate from the entries of $G^{-1}(x)$ (see Eq. (4.6)),

[TABLE]

Therefore eqs. (7.5) and (7.6) completely specify the inhomogeneous solution at order zero once the boundary constants $c_{i}^{(0)}$ are fixed, for instance by imposing the regularity of the solutions at specific kinematic points.

We have already observed that the system (3.7), or equivalently the third order differential equation (5.1), has regular singular points at $x=1$ and $x=\pm\infty$ , which correspond, respectively, to the pseudo-thresholds $s=4m^{2}$ and $s=0$ of the equal-mass banana graph. One can show that demanding the regularity of $\mathcal{I}_{i}^{(0)}(x)$ at such points is indeed sufficient in order to fix the three integration constants $c_{i}^{(0)}$ . In fact, by imposing regularity directly on the system of differential equations, one can determine three independent linear relations, which must be satisfied by the master integrals on the pseudo-thresholds. In particular, regularity at $x\to 1^{\pm}$ requires

[TABLE]

whereas at $x\to\pm\infty$ we find

[TABLE]

It is worths observing that, since $x\to\pm\infty$ corresponds to $s\to 0^{\pm}$ , the two conditions (7.8, 7.9) consistently reproduce the IBPs identities between the three-loop vacuum diagrams to which the master integrals are reduced in the zero-momentum limit.

It is particularly convenient to fix explicitly the boundary constants by working in the region $1<x<\infty$ , since the end-points of this region corresponds exactly to the two pseudo-threshold where we impose the regularity conditions (7.7-7.9). If we specify Eq. (7.5) to the interval $(1,\infty)$ and apply the rotation (7.6), we get

[TABLE]

where we have chosen as integration-base point $x_{0}=1$ and re-introduced the superscript $(1,\infty)$ for all quantities that require analytic continuation. We remark that, when applied to Eq. (7.10), the definition (7.6) of the function $\mathcal{R}^{(1,\infty)}_{i}(x)$ must be interpreted in terms of the homogeneous solutions $G^{(1,\infty)}(x)$ , which are defined in (A.4). Due to the choice of the integration base-point, in the $x\to 1^{+}$ limit all integrals appearing in the r.h.s of (7.10) vanish and the master integrals become

[TABLE]

The limiting behaviours of the homogenous solutions $H_{i}(x)$ , $J_{i}(x)$ , $I_{i}(x)$ at the two pseudo-thresholds are discussed in Appendix A and it is easy to verify that, when the expansions at $x\to 1^{+}$ (LABEL:eq:exp1p) are plugged into Eq. (7.11), the regularity constraint (7.7) is satisfied by demanding

[TABLE]

In a similar way, we can impose regularity at $x\to+\infty$ by making use of the expansions (A.30). Remarkably, due to the presence of logarithmic divergences $\ln(1/x)$ in the expansion of the homogenous solutions which must cancel in the expression of the master integrals, Eq. (7.8) allows us to fix at once $c_{1}^{(0)}$ and $c_{2}^{(0)}$ , which are given by

[TABLE]

In the second equalities we have simply performed the change of variable $t\to 1/y$ in order to map the integration range to $0<y<1$ . As a consistency check, we have verified that these values of the integration constants are consistent also with the regularity condition for the third master, given by Eq. (7.9).

Although we were not able to determine an analytic expression of the boundary constants, their representation as definite integrals (7.13) allows a high-precision numerical evaluation. We get for example

[TABLE]

The representation (7.10) of the master integrals, which is now fully determined, is valid for $1<x<\infty$ . The expression of $\mathcal{I}_{i}^{(0)}$ in the other kinematic regions (and in particular for $0<x<1/4$ , $s>16m^{2}$ , where the master integrals develop an imaginary part) can be obtained by analytic continuation of Eq. (7.10). The details of the analytic continuation are presented in Appendix A. As a summary of the results there discussed, let us just stress that the determination of a set explicitly real solutions $G^{(a,b)}(x)$ in each region and the study of their leading behaviour in the proximity of the singular points can be used to define, for any value of $x$ , a representation of the solutions which involves individually real-valued integrals only and, therefore, allows a fast and accurate numerical evaluation. A plot of the numerical results obtained through our representation compared against the computer code SecDec 3 [54] is shown in Fig. 3.

8 Conclusions

In this paper, we dealt with two important issues for the calculation of complicated multiloop Feynman integrals. The first, which has recently received a lot of attention, has to do with the explicit construction of a full set of homogeneous solutions of the system of couple differential equations satisfied by a family of Feynman integrals. The second, involves the classes of functions that are required for their calculation, beyond the well understood multiple polylogarithms. To address both issues in a non-trivial environment, we considered the calculation of the three-loop massive banana graph, which is the simplest graph known to satisfy an irreducible third-order differential equation. Complementing the findings reported in [40], we showed that all independent homogeneous solutions can be found by evaluating the maximal cut of the graph on the independent contours which do not cross any branch cuts of the integrand. We showed explicitly that in the case of the three-loop massive banana graph, there are only three possible choices of independent contours and that they indeed provide all three independent homogeneous solutions of its third-order differential equation. Our findings are in perfect agreement with the ones recently presented in [51].

Given the three independent homogeneous solutions as integrals over contours, we showed how to evaluate them explicitly in terms of products of complete elliptic integrals. The result, which was found long ago with very different methods by Joyce [49], is very interesting for two reasons. On the one hand, it is the first known generalization of the two-by-two differential equations satisfied, for example, by the two-loop sunrise graph, and which give rise to integrals over elliptic integrals. On the other hand, while generalizing the case above, it does not require to introduce truly new functions, as the homogeneous solutions can be still written as products of elliptic integrals only. While this is known to be an accident of the equal-mass case, it is nevertheless a very interesting laboratory to study the functions required to evaluate Feynman integrals which fulfil differential equations of order higher than two.

The results described in this paper show that the maximal cut can be used to solve complicated homogeneous differential equations, even when the number of coupled equations is higher than two. At variance with the $2\times 2$ case, which seems to invariably require complete elliptic integrals, $3\times 3$ and higher cases might require the introduction and classification of more complicated functions. We look forward to applying this method to increasingly complicated cases of relevance for high-energy physics phenomenology.

Acknowledgements

One of the authors (L.T.) is grateful to the Mainz Institute for Theoretical Physics (MITP), where some of the ideas at the basis of this project were developed. We wish to thank in particular Johannes Henn and Stefan Müller-Stach, whose encouragement and suggestions were crucial at different stages of the project. We also acknowledge many clarifying discussions with Claude Duhr, Harald Ita, Pierpaolo Mastrolia and Ettore Remiddi. Finally, we wish to thank David Broadhurst for interesting comments and for pointing us to the correct literature on cubic lattice Green functions. A.P. wishes to thank the Institute for Theoretical Particle Physics of the Karlsruhe Institute of Technology for hospitality during the development of this project.

Appendix A Analytic continuation

In this appendix we describe the analytic continuation of the master integrals derived in Section 7 to arbitrary values of $x=4m^{2}/s$ . We start from the analytic continuation of the homogenous solutions by first defining, for each kinematic region $a<x<b$ , a set of real-valued homogeneous solutions $G^{(a,b)}(x)$ and then by matching their limiting behaviours in order to link them across the singularities of the differential equation. Finally we will apply these results to Eq. (7.10) and obtain the analytic continuation of the inhomogeneous solution.

A.1 Homogenous solutions

In Sections 5-6 we have obtained, through different approaches, two different representation of the homogenous solutions, corresponding to Eq.s (5.46) and (6.14). Although the two representations have been shown to be completely equivalent, we decide to work with the latter, since it leads to more compact expression. Therefore, we consider homogenous solutions written in terms of products of complete elliptic integrals of with arguments

[TABLE]

The solutions (6.14) are explicitly real in the region $(1,\infty)$ . In order to define a set of real solutions in the other three regions, $(-\infty,0)$ , $(0,1/4)$ and $(1/4,1)$ , we can make use of well-known identities between complete elliptic integrals such as

[TABLE]

which establish linear relations between elliptic integrals with different reality domains. Therefore, by extending the set of building-blocks of the homogenous solutions to the elliptic integrals

[TABLE]

one can easily obtain, for each region $(a,b)$ , a matrix of homogeneous solutions $G^{(a,b)}(x)$ with real entries. In the following we list, for each region, one possible choice of real solutions for the first master integral, which correspond to the first row of $G^{(a,b)}(x)$ . As we have already observed, the other two rows can be obtained by applying the differential operators (5.2, 5.3) to the first one.

•

$-\infty<x<0$ :

[TABLE]

with Wronskian

[TABLE]

•

$0<x<1/4$ :

[TABLE]

with Wronskian

[TABLE]

•

$1/4<x<1$ :

[TABLE]

with Wronskian

[TABLE]

•

$1<x<\infty$ :

[TABLE]

with Wronskian

[TABLE]

A.2 Analytic continuation of the homogenous solution

Once real homogenous solutions $G^{(a,b)}(x)$ have been found in each region $(a,b)$ , we must match them across the four singular points $x=0,1/4,1$ and $x=\pm\infty$ in order to analytically continue the homogenous solutions the whole range $-\infty<x<\infty$ . Given the matrices $G^{(a,b)}(x)$ and $G^{(b,c)}(x)$ of real solutions defined in the adjacent intervals $(a,b)$ and $(b,c$ ), the analytic continuation amounts to define a matching matrix $M^{(b)}$ ,

[TABLE]

which allows to continue $G^{(a,b)}(x)$ to the region $b<x<c$ . The matrix $M^{(b)}$ can be obtained by assigning a small imaginary part to $x\to x-i0^{+}$ (the sign of which is inherited from the Feynman prescription $s\to s+i0^{+}$ ) and by equating the $x\to b^{+}$ limit of the two sides of (A.12).

This procedure leads the four matching matrices

[TABLE]

which, consistently with Eq. (A.12), satisfy

[TABLE]

The limits of the homogenous solutions (A.4 -A.10) close to the singular points, which have been used to obtain (A.26), can be easily calculated with the help of computer algebra system such as Mathematica and, therefore, we will not write them down explicitly. As an example, we will just list below the leading behaviour of the homogenous solutions (A.10) at the end-points of the region $(1,\infty)$ , which have been also used in Section 7 in order to fix the boundary constants of the inhomogeneous solution.

The limit of $G^{(1,\infty)}(x)$ for $x\to 1^{+}$ is

[TABLE]

where we have defined

[TABLE]

The leading behaviour of $G^{(1,\infty)}(x)$ for $x\to+\infty$ is instead

[TABLE]

A.3 Analytic continuation of the inhomogeneous solution

We are finally in the position to analytically continue the inhomogeneous solutions to arbitrary values of $x$ . As it is explicitly shown by Eq. (7.10), the inhomogeneous solution is defined, for $x>1$ , in terms of integrals of the functions $\mathcal{R}^{(1,\infty)}_{i}(x)$ which, in turn, depend on the homogenous solutions $G^{(1,\infty)}(x)$ and on their Wronskian $W^{(1,\infty)}(x)$ . Therefore, in order to extend the integral representation (7.10) to the other kinematic regions, it is sufficient to analytically continue the elements of $G^{(1,\infty)}(x)$ appearing in the definition (7.6) of $\mathcal{R}^{(1,\infty)}_{i}(x)$ , by making use of the matching matrices $M^{(b)}$ , as prescribed by Eq. (A.12). In this way all imaginary parts (whenever they are present) are made explicit and, as a result, we obtain a representation of the solution which involves, for any $x$ , the evaluation real integrals only.

We start by considering the analytic continuation to $1/4<x<1$ . The matching matrix $M^{(1)}$ , which has been defined in Eq. (A.26), can be used in order to express the homogeneous solutions $G^{(1,\infty)}$ in terms of the set of real solutions defined in $(1/4,1)$ ,

[TABLE]

In addition, it is easy to see that, with the Feynman prescription $x\to x-i\epsilon$ , the Wronskian is analytically continued for $1/4<x<1$ as

[TABLE]

By acting with Eq.s (A.31) and (A.32) on the inhomogeneous solution (A.26), we can write the master integrals in region $(1/4,1)$ in terms of individually real-valued integrals as

[TABLE]

where $\mathcal{R}^{(1/4,1)}_{i}(x)$ are combinations of homogenous solutions defined by

[TABLE]

We can now continue the solution to $0<x<1/4$ , where the master integrals develop an imaginary part. The region $(0,1/4)$ must be linked to $(1,\infty)$ by passing through region $(1/4,1)$ . This means that, according to the definition (A.12), the homogenous solutions $G^{(1,\infty)}(x)$ are continued in terms of the real-valued solutions defined for $0<x<1/4$ as

[TABLE]

where $M^{(1/4)}$ and $M^{(1)}$ are the matching matrices given in Eq. (A.26). In this case, Wronskian is trivially continued,

[TABLE]

and by making use of Eq.s (A.35) and (A.36) we can write master integrals in region $(0,1/4)$ as

[TABLE]

where the integration constants $b_{i}^{(0)}$ are defined by

[TABLE]

and $\mathcal{R}^{(0,1/4)}_{i}(x)$ are the combinations of homogenous solutions

[TABLE]

Finally, the expression of the master integrals in the Euclidean region $x<0$ can be obtained by matching the homogenous solutions at infinity, according to Eq. (A.12),

[TABLE]

with the matching matrix $M^{(\infty)}$ defined in Eq.(A.26). The Wronskian can be directly continued to negative values of $x$ ,

[TABLE]

and by acting again with (A.40) and (A.41) on Eq. (7.10), we obtain the expression of the master integrals in region $(-\infty,0)$ ,

[TABLE]

where $d_{1}^{(0)}$ is the integration constant

[TABLE]

and $\mathcal{R}^{(-\infty,0)}_{i}(x)$ are the combinations of homogenous solutions

[TABLE]

We stress here that similar results in the region $(-\infty,0)$ could be obtained by matching the solutions in $x=0$ . This point, nevertheless, is more delicate, due to the divergence in the Wronskian, i.e. $1/W(x)\sim 1/x^{3}$ as $x\to 0$ , and we preferred for this reason to continue passing through $x=\pm\infty$ .

Appendix B Proof of Eq.s (5.43)-(5.44)

In this appendix we give a brief derivation of Eq.s (5.43)-(5.44). We first reproduce the proof of Eq. (5.44), which was first presented in [53] (see Eq.(33) therein), and then we use similar arguments to derive Eq. (5.43).

The evaluation of Eq. (5.42) requires the study of the integral

[TABLE]

The analytic expression of (B.1) can be obtained by first studying the following auxiliary integral

[TABLE]

where $K_{0}(x)$ is the modified Bessel function of the second kind,

[TABLE]

which satisfies the identity

[TABLE]

Eq. (B.4) allows to trivially perform the integration over $z_{2}$ ,

[TABLE]

The integral over $z_{1}$ is now in standard form (see for instance Eq.2.16.36.2 of [55]) and it can be evaluated in terms of an elliptic integral of the first kind,

[TABLE]

from which we immediately see that Eq. (B.4) corresponds to the value of the auxiliary integral (B.2) at $\omega=0$ . In order to evaluate $I_{1}(0)$ , we go back to Eq. (B.2) we start by performing the $dt$ integration , for which we can use the distribution identity

[TABLE]

The term proportional to $\delta(z_{1}+z_{2})$ in the right-hand-side of Eq. (B.7) has no support in the region where $I_{1}(\omega)$ is defined, therefore we have

[TABLE]

which, if we set $\omega=0$ , reduces to

[TABLE]

This last integral is connected to the master formula (see Eq(3.3) of [56])

[TABLE]

where $I_{\mu}(z)$ is the modified Bessel function of the first kind, the function $W_{\mu}(k)$ is related to associated Legendre polynomial $P_{-1/2}^{\mu/2}$ ,

[TABLE]

and the arguments $k_{\pm}$ are defined by

[TABLE]

The expansion of Eq.(B.10) around $\mu=0$ allow us to express a set of integrals of three Bessel functions as a product of two complete elliptic integrals. In fact, by making use of

[TABLE]

one can easily check that that Eq. (B.10) implies

[TABLE]

Thanks to Eq. (B.15) we can finally evaluate $I_{1}(0)$ ,

[TABLE]

which proves Eq. (5.44).

The proof of Eq. (5.44), which requires the evaluation of the integral

[TABLE]

proceeds along the same lines. We start from the auxiliary integral

[TABLE]

which, by using in order (B.4), becomes

[TABLE]

As in the previous case, the integral over $z_{1}$ can be evaluated in terms of an elliptic integral of the first kind (see for instance Eq.2.16.36.2 of [55]),

[TABLE]

from which we see that Eq. (B.17) corresponds to $I_{2}(0)$ . Therefore, in order to determine the value of the auxiliary integral at zero, we first make use of (B.7) in Eq. (B.18) in order to integrate over $t$

[TABLE]

and then, after setting $\omega=0$ ,

[TABLE]

we make use Eq. (B.14) and obtain

[TABLE]

which proves Eq. (5.43).

Bibliography56

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] F. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions , Phys.Lett. B 100 (1981) 65–68.
2[2] K. Chetyrkin and F. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops , Nucl.Phys. B 192 (1981) 159–204.
3[3] A. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation , Phys.Lett. B 254 (1991) 158–164.
4[4] Z. Bern, L. J. Dixon, and D. A. Kosower, Dimensionally regulated pentagon integrals , Nucl. Phys. B 412 (1994) 751–816, [ hep-ph/9306240 ].
5[5] E. Remiddi, Differential equations for Feynman graph amplitudes , Nuovo Cim. A 110 (1997) 1435–1452, [ hep-th/9711188 ].
6[6] T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions , Nucl.Phys. B 580 (2000) 485–518, [ hep-ph/9912329 ].
7[7] J. M. Henn, Multiloop integrals in dimensional regularization made simple , Phys.Rev.Lett. 110 (2013) 251601, [ ar Xiv:1304.1806 ].
8[8] A. B. Goncharov, Geometry of configurations, polylogarithms, and motivic cohomology , Adv. Math. 114 (1995), no. 2 197–318.