Differential equations and dispersion relations for Feynman amplitudes. The two-loop massive sunrise and the kite integral

It is shown that the study of the imaginary part and of the corresponding dispersion relations of Feynman graph amplitudes within the differential equations method can provide a powerful tool for the solution of the equations, especially in the massive case. The main features of the approach are illustrated by discussing the simple cases of the 1-loop self-mass and of a particular vertex amplitude, and then used for the evaluation of the two-loop massive sunrise and the QED kite graph (the problem studied by Sabry in 1962), up to first order in the (d-4) expansion.


Introduction
In the last years we have assisted to an impressive increase in our knowledge of the mathematical structures that appear in multiloop Feynman integrals, thanks to the combined use of various computational techniques, such as to the method of differential equations [1][2][3], the introduction of a class of special functions (dubbed originally harmonic polylogarithms, HPLs [4,5], they came out to be a subset of the much larger class of multiple polylogarithms, MPLs, see [6][7][8][9] and references therein), the definition of a so-called canonical basis [10] for dealing with increasingly larger systems of differential equations and the use of the Magnus exponentiation [11].
However, most of the above results have been obtained in the massless limit; indeed, the situation for massive amplitudes is different, as the two-loop massive sunrise (which has three propagators only) is still the object of thorough investigation [12][13][14][15][16][17][18][19][20]. A general approach to the study of arbitrarily complicated systems of differential equations within difference field theory has been recently proposed in [21].
In this paper we will show that the study of the imaginary parts and related dispersion relations satisfied by the Feynman amplitudes, within the differential equation frame, can provide another useful practical tool for their evaluation in the massive case as well.
The imaginary parts of Feynman graphs can be obtained in various ways. To start with, one can use Cutkosky-Veltman rule [22][23][24] for integrating directly the loop momenta in the very definition of the graphs. When the d-continuous dimensional regularization is used, nevertheless, that is practical only in the simplest cases. Another possibility is the extraction of the imaginary part from the solution of the differential equations, which of course requires the knowledge of the solution itself. More interestingly, one can observe that often the differential equations become substantially simpler when restricted to the imaginary part only, so that their solution can become easier.
In any case, once the imaginary part of some amplitude A(d; u), say ImA(d; u), is obtained, one has at disposal the dispersive representation for A(d; u), namely an expression of the form (where the limits of integration have been skipped for ease of typing). Such a representation turns out to be very useful when the amplitude A(d; u) appears within the inhomogeneous terms of some other differential equation, regardless of the actual analytical expression of A(d; u). Indeed, as the whole dependence on u is in the denominator (t − u) one can work out its contribution by considering only that denominator, freezing, so to say, the t -integration and the weight ImA(d; t) until the dependence on the variable u (the variable of the differential equation) has been properly processed. Let us emphasize, again, that such a processing is, obviously, fully independent of the actual form of ImA(d; t).
In the following, we will illustrate the above remarks in a couple of elementary applications and then use them in the case of the two-loop QED-kite, i.e. the two-loop electron self-mass in QED, already studied by Sabry [25] long ago. The study of the kite amplitudes requires in turn the knowledge of the two-loop massive sunrise, which appears as inhomogeneous terms in their differential equations. Indeed, the imaginary part [26] and related dispersion relations [27,28] have been already exploited long ago for studying the zeroth order of the sunrise and the kite integral. In this paper our goal is more general, as we will show how to use them consistently within the differential equations approach, which will allow us to investigate the solution at any order in the (d − 4) expansion.
The paper is organized as follows. We begin in section 2 studying the imaginary part of the one-loop self mass and its dispersion relation for generic values of the dimensions d. We elaborate on its calculation both from Cutkosky-Veltman rule and from the differential equations. In section 3 we study a particular vertex amplitude through the differential equations method. The one-loop self-mass appears as inhomogeneous term in the equations and we show that their evaluation can be simplified, once the one-loop self-mass is inserted as dispersive relation. In section 4 we move our attention to the two-loop sunrise graph, which we write as iteration of two one-loop bubbles. This allows us to derive an extremely compact representation valid for generic d, from which one can show that, at every order in (d − 2) (and therefore also in (d − 4)), the sunrise can be written as a one-dimensional integral over a square root of a quartic polynomial, times a combination of multiple polylogarithms only. The simplicity of this result motivates us to look more systematically for a similarly simple structure using differential equations from the very beginning for the whole integral family of the kite. In section 5 we discuss the notation and describe the master integrals which have to be computed. In section 6 we provide the solution of the simple topologies, which can be written in terms of HPLs only. Then in section 7 we start a systematic study of the differential equations of the sunrise graph. It is known that the solution for the sunrise graph is somewhat simpler when its Laurent series is considered in (d − 2) instead of in (d − 4); however, we find more convenient to expand all the master integrals in (d − 4) from the very beginning. To that aim, by using the well known fact that any Feynman integral in d − 2 dimensions can be written as a linear combination of integrals in d dimensions, we build up a new, equivalent basis of master integrals for the sunrise whose expansion in (d − 4) is identical to the expansion of the original masters in (d − 2). Once we have a convenient basis and the corresponding differential equations, we show how to solve them iteratively in section 8. We conclude the section providing explicit analytical results for both master integrals for the first two non-zero orders and showing how to extract their imaginary parts and write dispersion relations for them. We move then to the kite integral in section 9, where we show how the representation of the sunrise as a dispersion relation is particularly convenient, as it allows to write a compact solution for the first two orders of the kite integral. Finally we conclude in section 10. We enclose different appendices where we provide further mathematical details and explicit derivations.

The 1-loop self-mass: imaginary part and dispersion relation
We define the integration over a loop momentum k in d continuous dimensions as so that the tadpole amplitude Tad(d; m) reads .
We then consider the 1-loop "bubble" . (2.4) We work in the Euclidean metric such that q 2 is positive when q is spacelike. At q = 0 one has at once Cutkosky-Veltman rule gives for the imaginary part of the bubble amplitude in d-continuous dimensions and for s = −q 2 > (m 1 + m 2 ) 2 the expression where we introduced the usual Källen function and the d dependent coefficient whose expansion for d ≈ 2 reads As a further remark, Eq. (2.6) can also be written as Once the imaginary part is given, we can write a dispersion relation for the one-loop bubble where Bub(d; 0, m 1 , m 2 ), which is given in Eq. (2.5), contains a pole at d = 4, while the integral is convergent for d < 6. The 1-loop self-mass amplitude Eq. (2.4), which in the following will be written as Bub(d; s) for ease of typing, is known to satisfy the following differential equation in s (2.13) where the inhomogeneous term, N(d; s) is given by (2.14) The homogeneous equation associated to Eq. (2.13) is (obviously) One sees immediately that ImBub(d; s), Eq. (2.10) satisfies the homogeneous equation, for any value of d. That fact is hardly surprising, yet it deserves some comments. When looking for a solution of an equation like Eq. (2.13), it can be convenient, in order to fix the boundary conditions, to start by considering values of the variable s for which the solution is expected to be real (typically, s = 0 or s negative, i.e. in the spacelike region). But as one is also interested in the value of the solution for timelike, physical values of s, one is naturally lead to consider the solution as a complex analytical function of the argument s, to be evaluated along the whole line s + i , with s real and varying in the range −∞ < s < +∞ and small and positive (the Feynman prescription). As the singular points of the equation correspond to real values of s, such as for instance s = (m 1 ± m 2 ) 2 , the function has no singularities along the s + i line, so that its value is fully determined by the analytic continuation in terms of the initial boundary conditions. Moreover, one might be interested in considering separately the real part ReBub(d; s) and the imaginary part ImBub(d; s) of the solution Bub(d; s) = ReBub(d; s) + iImBub(d; s). In so doing, as the inhomogeneous term is real, Eq. (2.13) splits into the two equations The evaluation of the solutions in the various regions is almost trivial, but one needs the knowledge of three constants, c 1 , c 2 , c 3 to actually recover the imaginary part Eq. (2.10) (in that case the constants are, obviously, c 1 = c 2 = 0, c 3 = πB d /2). Summarizing, the evaluation of the imaginary parts alone within the differential equation approach is much simpler than the evaluation of the complete solution (real and imaginary parts), but requires some additional external information (such as the knowledge of the regions in which the imaginary part vanishes and its normalization when not vanishing).

The 1-loop self-mass and the 1-loop equal mass triangle
In this section we show how to use Eq. (2.11) in the solution of the differential equation for a particular massive triangle amplitude, namely with q = p 1 +p 2 , p 2 1 = p 2 2 = 0 and −q 2 = s (q 2 is positive when q is spacelike). The differential equation in the variable s for the amplitude Tri(d; s) reads We can then use Euler's method to write the general solution for the triangle as follows where c(d, m) is an integration constant, depending in general on d and m. We can fix the integration constant requiring that for s → 0 the amplitude is not divergent, which implies The tadpole, Eq. (2.3), is of course independent of u and for the bubble we use its dispersive representation Eq. (2.11) Let us recall that the above integral is convergent for d ≈ 2 and that if one is interested in d ≈ 4, one can use its subtracted version Eq. (2.12). Assuming for definiteness s < 4m 2 , and therefore ignoring for the moment the +i prescription in Eq. (3.6), the triangle amplitude becomes The integration in u is trivial and we get Note that the above result holds for any d (within the considered range) independently of the actual explicit form of the inserted amplitude Bub(d; u). If 0 < s < 4m 2 the result (3.8) is real, while for s > 4m 2 it develops an imaginary part. In order to properly extract it, it is enough to notice that, for s > 4m 2 , the s → s + i prescription gives and ∞ 4m 2 (3.10) Collecting results and combining the various terms, the imaginary part of Tri(d; s) for s > 4m 2 becomes (3.11) It is to be noted, again, that the above result has been obtained from Eq. (3.8) independently of the explicit analytic expression of ImBub(d; u). As a check, we can write the dispersion relation for the triangle amplitude in terms of its imaginary part (we take s < 4m 2 for simplicity) By exchanging the order of integrations according to Summarizing, the use of the dispersive representation of the inserted amplitude Bub(d; u) in the Euler form of the solution of the differential equation for the triangle amplitude gives, almost at once, the explicit form of Tri(d; s), Eq. (3.8) in terms of ImBub(d; u). The imaginary part ImTri(d; s) of the triangle amplitude Eq. (3.11) can also be written in terms of ImBub(d; u), without explicit reference to the analytic form of the latter. The resulting dispersion relation Eq. (3.12) can be useful if Tri(d; u) appears within the inhomogeneous terms of the equations for the amplitudes of some other process (such as for instance the QED light-light graphs).

The sunrise as iteration of the bubble graph
Let us consider the sunrise scalar amplitude defined as , (4.1) where the integration measure is defined in Eq. (2.1) and we work in the Euclidean metric for simplicity. It is well known that the sunrise graph with different masses possesses four master integrals, which reduce to two in the case of equal masses [12]. In this section we will not try to give a full solution for all the masters integrals, but instead we will limit ourselves to considering the scalar integral (4.1) only and try to study its iterative structure in (d − 2), or equivalently in (d − 4). One possible way to do this is by noting that the sunrise integral can be written as , (4.2) and according to Eq. (2.4) the integral in the momentum l is simply a one-loop bubble with masses m 2 and m 3 and momentum q = (p − k), .
The dispersive representation Eq. (2.11) then gives with ImBub(d; t, m 2 , m 3 ) given by Eq. (2.6). As q = p − k, Eq. (4.2) becomes . (4.5) Now clearly the integral in k can be seen again as a one-loop bubble amplitudes, this time with (squared) masses m 2 1 and t . Using again the formula (2.11) we get We can obtain an equivalent representation by further exchanging the integrations in the variables t and v such that, by recalling Eq. (2.6), we are left with (4.7) Eq. (4.7) is the main result of this section. From it we obtain at once, Note that Eq. (4.8) is nothing but the d-dimensional three-body massive phase space and Eq. (4.7) could indeed have been obtained also by computing first the imaginary part of the sunrise graph using Cutkosky-Veltman rule, and then writing a dispersion relation for it. Remarkably, the complexity of the result in the general mass case is practically the same as in the equal mass case Let us further emphasize that Eqs. (4.6), (4.7) and (4.8) are all true for generic, continuous values of d. Furthermore, their expansion in (d − n), where n is virtually any positive integer (and in particular in (d − 2)), is completely straightforward and generates only products of logarithms. 1 This implies in turn that, at every order in (d − 2), the integral in v in Eq. (4.6) can always be performed in terms of multiple polylogarithms only. This shows that, at every order in (d − 2), the sunrise integral can be written as a one-fold integral over the root of a quartic polynomial, times combinations of multiple polylogaritms. The result is interesting and it resembles similar results found for the finite term of a completely unrelated massless double box in N = 4 [29,30]. 2 Finally, the relation of this representation of the imaginary part of the sunrise Eq. (4.8) with the results obtained by the explicit solution of the system of differential equations for the two amplitudes of the sunrise problem (which involves two pairs of solutions, i.e. four functions altogether, see for instance section 8 of this paper) is also intriguing, but will not be further investigated here. Starting from the next section we will instead focus on the more general problem of computing the full set of master integrals of the kite graph using the differential equations method.

The differential equations for the kite master integrals
Let us consider the family of the integrals of the QED kite graph with three massive propagators and two massless ones, defined as I(n 1 , n 2 , n 3 , n 4 , n 5 where dashed lines represent massless propagators. The five denominators are chosen as with −p 2 = s and p 2 > 0 when p is spacelike. The integration measure is defined as in Eq. (2.1) such that according to Eq. (2.3) the one-loop tadpole reads .
The integral family (5.2) can be very easily reduced to master integrals using, for example, Reduze 2 [31,32]. In order to simplify the notation we put m = 1 and define u = s/m 2 . We find 8 independent master integrals which we choose as follows Most of the master integrals are very simple and have been already studied thoroughly in the literature. In particular M 1 ,. . . ,M 5 are known and can be written in terms of HPLs only. The remaining three integrals, M 6 , M 7 and M 8 , cannot be expressed in terms of MPLs and will be the main topic of this paper. Note that M 6 and M 7 are the two master integrals of the two-loop massive sunrise with equal masses, see Eq. (4.1). As we will see, M 6 and M 7 satisfy a system of two coupled differential equations, with M 6 appearing further within the inhomogeneous terms of the differential equation for M 8 .
As usual, we are interested in the Laurent expansion of the master integrals for d ≈ 4. The computation of the first five integrals in terms of HPLs is straightforward. In particular, it can be simplified by the choice of a canonical basis in d ≈ 4, which can be found following the methods described in [33,34]. For the last three integrals, instead, a canonical basis in the usual sense cannot be found and we will have to resort to different arguments in order to put the system of differential equations in a form that is suitable for their integration. We choose the following canonical basis for the simple topologies while for the non-trivial topologies we introduce The system of differential equations for the first five masters integrals can then be written as where the matrix A(u) reads (5.8) The differential equations for the last three integrals cannot be put in a similarly simple form and we will write them explicitly later on, once we come to study them.

The simple kite master integrals
Let us focus on the first five integrals. If we start from the differential equations (5.7), carrying out the integration in terms of harmonic polylogarithms is straightforward. As it is well known, harmonic polylogarithms are a special case of multiple polylogarithms and for convenience of the reader, we recall here their iterative definition. We start at weight one defining The multiple polylogarithms are then iteratively defined at weight n as follows G(a 2 , . . . , a n ; x) . (6.2) Note that all integrals have a cut at s = m 2 , u = 1, i.e. they are real for u < 1 and develop an imaginary part for u > 1 whose sign is fixed by Feynman's prescription u → u + i 0 + . We present here the solution valid for 0 < u < 1. The analytic continuation to the physical region can be then easily obtained by continuing to u > 1 with u → u + i 0 + . Thanks to the choice of a canonical basis the solution takes a particularly compact form and all sub-topologies, up to weight 4, fit in one single page. (6.7)

The choice of the basis for the sunrise amplitudes
We move now to consider the last three integrals. First of all we need to focus on the two master integrals of the two-loop sunrise graph, i.e. f 6 (d; u) and f 7 (d; u). They satisfy a system of two coupled differential equations Let us recall here that amplitude f 1 (d; u) appearing within the inhomogeneous term corresponds to the product of two tadpoles and is in fact constant, according to Eq. (6.3).
Using the methods described in [35] one can show that it is not possible to decouple the system, in any even number of dimensions d = 2 n, n ∈ N, by taking simple linear combinations of the masters integrals with rational coefficients. Indeed, it is very well known that the solution cannot be expressed in terms of MPLs only and elliptic generalizations of the latter must be introduced [13][14][15][16][17][18].
In order to simplify the integration of these two integrals we will proceed as follows. We will start considering the integrals for d ≈ 2. The reason for this is two-fold. First, when d = 2 the two master integrals f 6 (2; u) and f 7 (2; u) are finite. Second, as we will see explicitly, their imaginary parts when d = 2 are particularly simple. That is important because the imaginary parts are related to the solutions of the corresponding homogeneous system, which in turn are the building blocks for the iterative solution of the 2 × 2 differential system (7.1) through Euler's method. Those considerations will allow us to determine a basis of master integrals for which we can easily solve the differential equations as a Laurent series in (d − 2). At this point, we could solve the system as Laurent series in (d − 2) and then use Tarasov shifting identities [36] in order to obtain the corresponding coefficients of their Laurent series in (d − 4), which are the physically relevant results. Instead of proceeding in this way, though, we will use the technique described in [35] (see appendix B therein) in order to build up a new basis of master integrals which fulfills the very same differential equations, but this time with d → d − 2. This implies that the series expansion in (d − 4) of the new basis will be formally identical to that of the former basis in (d − 2). This will allow us to treat more consistently everything in d ≈ 4 from the very beginning.
7.1. Simplifying the differential equations in d = 2 The system of differential equations (7.1) has four regular singular points, i.e. u = 0, u = 1, u = 9 and u = ±∞, we will therefore need to consider the solution in the four different regions Physically, the point u = 9, s = 9 m 2 , corresponds to the three massive particle cut and we expect the master integrals to develop an imaginary part as u > 9. Now, since the tadpole does not have any cut in u, the imaginary parts of the master integrals f 6 (d; u) and f 7 (d; u) must satisfy the associated homogeneous system. We have already computed the imaginary part of the first master integral in section 4 for generic d. For d = 2, a straightforward application of Cutkosky-Veltman's rule to f 7 (2; u) as well gives 1 π Imf 6 (2; u) = I (0, u) where the functions I (n, u) are defined as and R 4 (d, u) is the fourth-order polynomial Some of the properties of these functions are discussed in Appendices A and D. In particular, there it is shown that all functions I (n, u) can be expressed in terms of two independent functions only, say I (0, u) and I (2, u), which can be in turn expressed in terms of the complete elliptic integrals of first and second kind, and can be therefore considered as known analytically. Eqs. (7.2) suggest to perform the following change of basis so that obviously the imaginary parts of the functions g 6 (d; u) and g 7 (d; u) in d = 2 become 1 π Im g 6 (2; u) = I (0, u) , 1 π Im g 7 (2; u) = I (2, u) . (7.6) Note that these relations are true in d = 2 and do not change if we modify (7.5) by a term proportional to (d − 2). This freedom can be used to get rid of the term proportional to (d − 2) 2 in the second of Eqs. (7.1). While this is not strictly required, it indeed helps in simplifying the structure of the solution. We modify therefore Eq. (7.5) as follows where C(u) is a function of u only, to be determined by imposing that the term proportional to (d − 2) 2 in Eqs. (7.1) is zero. By writing down explicitly the differential equations for g 6 (d; u) and g 7 (d; u) one easily finds that there are two values of C(u) which would eliminate the unwanted term, namely At this level, there is no reason to prefer one choice over the other, we choose therefore C(u) = 6(u − 1), since this produces the most compact results. With this choice the differential equations become d du where we used f 1 (d; u) = 1, Eq. (6.3). The system can be written in matrix form as follows d du where the two matrices B(u), D(u) are defined as . (7.11) In order to be able to solve (7.9) as a Laurent series in (d − 2), as a first step we need to solve the homogeneous system for d = 2, i.e. we need to find a pair of two solutions, say (I 1 (u), I 2 (u)) and (J 1 (u), J 2 (u)), such that the matrix of the solutions Note in particular that since Tr(B(u)) = 0 , (7.14) the Wronskian of the four solutions, W (u) = I 1 (u) J 2 (u) − I 2 (u) J 1 (u), must be independent of u. From its very definition, we find = Tr(B(u)) det (G(u)) = 0, (7.16) and W (u) must be a constant. This property is of fundamental importance to simplify the iterative solution of the differential equations, as we will see later on.

The choice of the basis for the expansion in (d − 4)
In the previous section we showed how to choose a basis of master integrals for the sunrise graph, whose differential equations take a particularly convenient form as far as their Laurent series in (d − 2) are considered. Since, as it is well known, any Feynman integral in d − 2 dimensions can be expressed as a linear combination of Feynman integrals in d dimensions, by following the method described in appendix B of [35] we define a new basis of master integrals by shifting (7.7) 1 (d, u)) as follows It is straightforward by direct calculation, and using (7.9), to prove that the new basis (7.18) satisfies the new system of differential equations d du As expected the system (7.19) is identical to (7.9), upon the formal substitution d → d − 2. This also implies that all the properties fulfilled by g 6 (d; u) and g 7 (d, u) in the limit d → 2, are also fulfilled by h 6 (d; u) and h 7 (d; u) in the limit d → 4. In particular the new master integrals are finite in d = 4 and their imaginary parts read 1 π Im h 6 (4; u) = I (0, u) , 1 π Im h 7 (4; u) = I (2, u) , (7.20) as Eqs. (7.6).

The solution of the differential equations
In this section we will show how to build the complete solution for the sunrise master integrals up to any order in (d − 4). We will solve the system as Laurent series in (d − 4). In order to do this, we first need to find the homogeneous solution in the limit d → 4, such that we can then use Euler's method of variation of constants in order to build up the complete non-homogeneous solution.

The homogeneous solution
As explained above, as a first step we need now to find two independent pairs of solutions for the homogeneous system associated to (7.19) d du The discussion in previous section already suggests how to find the first of the two pairs. Taking the imaginary part of (7.19) at d = 4 gives at once d du , (8.5) as expected.
In order to find a second pair of solutions, we go back to the definition of the functions I (n, u), introduced in Eq. (7.3) as the definite integral in b of the square root of the fourth-order polynomial R 4 (b, u), Eq. (7.4), between two adjacent roots. Since R 4 (b, u) has 4 roots, we are naturally brought to consider two similar sets of functions, defined by integrating between the other pairs of adjacent roots, say More details on these functions are provided in Appendix A. In particular, one can show that, also in this case, there are two "master integrals" for each set of functions, say J (0, u), J (2, u) and K(0, u), K (2, u). Moreover, one can show that the functions K(n, u) are not independent from the functions J (n, u) and we can therefore neglect them. We pick for definiteness J (0, u) and J (2, u) and we compute their derivatives finding J (0, u) and J (2, u), as they stand, are not solutions of the homogeneous system; it is nevertheless very easy to use them in order to build a proper solution for the system. Consider the new functions defined as By using (8.7) it is trivial to verify that their derivatives read (8.8) so that we can define our second set of solutions, again for 9 < u < ∞, as Summarizing we have found two pairs of independent real valued solutions, valid in the range 9 < u < ∞, such that their matrix We can now proceed and study their limiting behaviour on the two boundaries, i.e. u → 9 + and u → +∞. For u → 9 + we find (keeping only the leading logarithmic behaviour) On the other hand for u → +∞ we find As stated previously, the Wronskian (7.15) must be independent of u; when can computing it using any of the limits above, we find and as expected.
The solution described here is valid above threshold, i.e. for u > 9, but it is straightforward to extend those results and build up a complete set of solutions valid in the remaining regions, i.e. −∞ < u < 0, 0 < u < 1 and 1 < u < 9. The details are worked out explicitly in Appendix B. We end up in this way with 4 different matrices of real solutions G (a,b) (u), each valid in the interval a < u < b, and which can be continued from one region to the other using the matching matrices given in the same appendix, see in particular Eqs. (B.16) and (B.17). Note that in (8.14) and (8.15) we computed the value of the Wronskian in the region 9 < u < ∞, but we can normalize the solutions in the remaining three regions such that the same remains true in every interval (a, b), see Appendix B, (8.16)

The non-homogeneous solution
Once we have the homogeneous solution of the system for d = 4 we can use Euler's method of the variation of constants in order to write the complete solution of the system Eq. (7.19). The manipulations performed here are the same for all the regions a < u < b, we will therefore drop the superscripts (a, b) from all formulas for simplicity, writing for instance G(u) instead of G (9,∞) (u) etc. It will be then simple to specialize the results to the region of interest by picking the suitable set of solutions G (a,b) (u), see for instance Eqs. (8.10), (8.11) for the notation. We perform the rotation Thanks to the condition on the Wronskian (8.16), inverting the matrix G(u) is straightforward and we get We write therefore the system as d du where we introduced the matrix . (8.20) Written in this form, the iterative structure of the solution in powers of (d − 4) becomes manifest. The entries of the matrix M(u) read they contains rational functions and products of pairs of homogeneous solutions, i.e. products of complete elliptic integrals. It should be recalled at this point that not all products are actually linearly independent; because of the condition on the Wronskian, in fact, only one of the two combinations I 1 (u)J 2 (u) or I 2 (u)J 1 (u) is really independent, while the other can be removed using (8.16). Moreover, notice that the functions I k (u) and J k (u) fulfill the same differential equations, i.e. (8.4) or (8.8). By inverting them one can, for example, get rid of I 2 (u) and J 2 (u) in favor of I 1 (u) and J 1 (u) and their derivatives Substituting these relations into (8.21) and rearranging the terms, the matrix can be written in a much more compact form as where we used , (8.25) which can be easily proved starting from the condition on the Wronskian (8.16). Equations (8.24) are particularly interesting, as they show that the matrix M(u) can be written as a total differential of simple logarithms plus three new functions which are given by products of complete elliptic integrals and a polynomial in u. Once appropriate boundary values are known, the integration of the system (8.19) as a Laurent series in (d − 4) becomes, at least in principle, straightforward in terms of iterated integrals with the kernels given by the entries of M(u) (8.24). Given this result, it is indeed very tempting to try and define a new generalized alphabet composed by the six generalized letters appearing in Eq. (8.24). We will resist the temptation for now, and instead go ahead and see how far we can get with what we have. We will work for simplicity in the region 0 < u < 1 and use everywhere the solutions of the homogeneous system valid in this region, G (0,1) (u), see Appendix B. Working in 0 < u < 1 is also very convenient since we can easily fix the boundary conditions imposing the regularity of the two original master integrals, h 6 (d; u) and h 7 (d; u), at u = 0 + and u = 1 − . From Eqs. (7.19) we can read off the two conditions Having determined (8.26), we can now proceed with the integration of the differential equations. We start from (8.19) and expand everything in (d − 4) as follows (8.27) such that at order zero the equations reduce to d du , (8.28) and at first order we have instead d du , (8.29) where the previous order appears as inhomogeneous term. Note that this structure remains true at every order n, with n ≥ 1 d du In the next two sections we describe the integration of (8.28) and (8.29), which will allow us to write a compact result for the master integrals of the two-loop massive sunrise graph, h 6 (d; u) and h 7 (d; u) up to first order in (d − 4).

The two-loop massive sunrise at order zero
The integration of the order zero, Eq. (8.28), can be carried out simply by quadrature. Specializing formulas above in the region 0 < u < 1 we find The constants c (0) 6 and c (0) 7 can be fixed imposing (8.26). Note that the latter must be imposed on the original master integrals h 6 (d; u) and h 7 (d; u) and not on m 6 (d; u) and m 7 (d; u), with the relation between the two sets of functions given by Eq. (8.17). Expanding also the original masters integrals as Using the limiting values given in Appendix B and the definite integrals of Appendix C we obtain where the last integral can be performed by standard techniques using the integral representation for J (0,1) 1 (u), see for example [13]. We recall here the definition of the Clausen function Finally putting everything together we get For convenience, we provide here the limiting values of the master integrals in the two matching points u = 0 + and u The solution (8.35) is valid for 0 < u < 1. We can use the matching matrices defined in Appendix B to continue the solution in any other region. In particular, it is interesting to study the continuation above threshold, i.e. for u > 9 where the master integrals develop an imaginary part. By straightforward use of the formulas in the appendix we find, for 9 < u < ∞, (u) , (8.38) such that, as expected, the imaginary parts of the two master integrals at order zero in (d − 4) read Note that the simplicity of the imaginary part above threshold, u > 9, and the absence of an imaginary part for the intermediate region 1 < u < 9 is true only for the physical masters integrals h 6 (d; u) and h 7 (d; u). The rotated functions, m 6 (d; u) and m 7 (d; u), have no direct physical meaning and cannot be expected in general to develop an imaginary part only above the u > 9 threshold.
Having the imaginary part, we can write an alternative representation for the solution (8.35) as a dispersion relation where for h (0) 7 (u) we have used a doubly subtracted dispersion relation and fixed the boundary terms matching (8.40) to (8.35) for u = 0 + and u = 1 − . As showed in section 3, this representation is also particularly convenient if we need to integrate once more over it, for example whenever the sunrise appears as subtopology in the differential equations of more complicated graphs, see section 9.

The two-loop massive sunrise at order one
The order zero of the sunrise graph is special since its inhomogeneous term is very simple. In order to understand the general structure we want to integrate Eq. (8.29), which implies integrating over the matrix M(u) in (8.21), using (8.31) as inhomogeneous term. Again specializing the formulas for 0 < u < 1 and integrating by quadrature we get m (1) At this point the solution is written as a double integral over known functions. The entries of the matrix M(u) are rather complicated, see (8.21). Nevertheless the result can be greatly simplified using integration by parts identities and the condition on the Wronskian (8.16). By direct inspection of the matrix (8.21) it is clear that, at any order (d − 4) n , the result can only contain at most the following integrals where F (t) is a generic function of t and, at order n, it contains the order (n − 1) of the Laurent expansion of the functions m 6 (d; u) and m 7 (d; u). For n = 1, which is the case we are interested in, inspection of Eqs. (8.41) shows that F (t) is either a constant, or it can be one of the two functions Using integration by parts identities together with the condition on the Wronskian (8.16), one can show that not all integrals (8.42) are independent. 3 In particular one can re-express all integrals containing the products I 1 (t) J 2 (t) and I 2 (t)J 1 (t) with the rational prefactors appearing in (8.42), in terms of the remaining integrals with I 1 (t)I 1 (t), I 1 (t)J 1 (t) and J 1 (t)J 1 (t) only. This allows to substantially simplify the resulting expressions and, notably, eliminate all occurrences of double integrals over the products of functions I k (t) and J k (t), at the price of introducing simple logarithms -a non-trivial result. For simplicity we provide here the analytical expressions for the physical master integrals only, i.e. h (1) 6 (u) and h (1) 7 (u), omitting the intermediate ones for m (1) 6 (u) and m (1) 7 (u). The latter can anyway easily be recovered by rotating the solution back through the matrix G −1 (u), see Eq. (8.17).
Again in the region 0 < u < 1 we can easily fix the boundary values using the results of Appendices B and C and we find 4 h (1) Of course, one also needs to make use of the differential equations satisfied by the I k (t) and J k (t) in order to re-express the derivatives dI k (t)/dt and dJ k (t)/dt in terms of the I k (t) and J k (t). 4 As discussed in Appendix C, we do not present explicitly all integrals required to fix all limits. The complete list of definite integrals can be obtained by the authors.
where, again, the dispersion relation for h (1) 7 (u) is doubly subtracted in u = 0. It is clear that, at least in principle, the techniques described here for the integration of the first two orders of the two-loop massive sunrise, can be used also for the higher orders. The formulas are of course more cumbersome and, in general, it is not granted that the result can always be written in terms of one-fold integrals only, as for order zero and one, like in Eqs. (8.35), (8.44) and (8.45). Nevertheless one can show that, by using integration by parts as we did for the order one, also the order (d − 4) 2 can be substantially simplified.
One last comment is in order. The basis of master integrals that we have been considering, h 6 (d; u) and h 7 (d; u), was build by the shift d → d − 2 of the previous basis g 6 (d; u), g 7 (d; u), see Eq. (7.7). That implies that if we expand the latter as Laurent series in (d − 2)

The solution for the kite integral
As a last step we will use the results of the previous sections in order to write compact expressions for the first two non-zero orders of the kite integral. We will do this using the method sketched in section 3, namely we will derive the differential equations for the kite integral and then we will insert into it the solution for the sunrise graph given as a dispersive relation, see Eqs. (8.40), (8.52) and (8.53). We start by writing the differential equations for the master integral f 8 (d; u), defined in (5.6), while for the sunrise we use the modified basis defined in (7.18). The differential equation reads (9.1) Two properties are worth noticing in Eq. (9.1). First, only one of the two master integrals of the sunrise subgraph appears, namely h 6 (d; u). Second, it appears multiplied by a factor (d − 4) 3 . This fact is a consequence of the normalization adopted in (5.6), where, attempting to build up a basis similar to a canonical one, we rescaled all master integrals of suitable powers of (d − 4).
Note however that, even if in (5.6) also f 6 (d; u) and f 7 (d; u) are rescaled by (d − 4) 2 , one should recall that the masters integrals that we are effectively calculating for the sunrise graph (and which enter in the differential equation for the kite) are not f 6 (d; u) and f 7 (d; u), but instead h 6 (d; u) and h 7 (d; u), as defined in (7.18). The latter are obtained shifting (5.6) from d → d − 2, such that the factor (d − 4) 2 in front of f 6 (d; u) and f 7 (d; u) becomes effectively a (d − 6) 2 . In order to make the equations more symmetric, we could have therefore rescaled also h 6 (d; u) and h 7 (d; u) by (d − 4) 2 , reabsorbing in this way the corresponding factor in (9.1). We preferred, nevertheless, not to do that in order to avoid the confusion of one more change of basis. With the present normalization, the sunrise integrals start at order zero in (d − 4), which shows that their first contribution to the Laurent expansion of the kite integral is at order (d − 4) 3 .
We can now move to the actual integration of the equations. Once more we work in the region 0 < u < 1, where the boundary condition can be read off directly from Eq. (9.1), imposing regularity of f 8 (d; u) on the pseudo-threshold u = 0. This condition implies It is easy to see that f 8 (d; u) is finite in d → 4 and therefore its Laurent expansion reads 1, u) .
At this point one could, in principle, plug in the solution for the sunrise integral as given by (8.38). That introduces anyway unneeded complications. The easiest way to proceed is instead to insert the dispersive solution, Eq. (8.40). Upon doing this, the integration in u becomes straightforward in terms of multiple polylogarithms and, after fixing the boundary condition, we are left (somewhat surprisingly!) with an extremely compact result The very same exercise can be repeated for the next order, making use of the dispersion relations derived for the sunrise graph at order one (8.52), and of the previous order just computed (9.7). By integrating the differential equation and fixing the boundary condition we get f (4) where l (t) is defined in (8.51). Note that Eq. (9.8) contains a combination of, on one side, polylogarithms of weight 4 and, on the other, of integrals over elliptic integrals and polylogarithms of weight 2.

The analytic continuation of the solution
Here we want to show that also in this case the analytic continuation of our solution, Eqs. (9.7) and (9.8), is completely straightforward in the whole range −∞ < u < +∞. The kite integral has a first cut at u = 1 corresponding to s = m 2 , where the harmonic polylogarithms develop an imaginary part. The second cut is at u = 9, s = 9 m 2 , and the elliptic integrals develop further imaginary parts which can be easily computed using the results of Appendix B. Let us consider for example the first non-zero order, Eq. (9.7), and let us continue it in the two physically relevant regions, i.e. for 1 < u < 9 and then above the three-mass threshold 9 < u < ∞.
In this region the HPLs develop an imaginary part, whose sign is fixed by Feynman's prescription u → u + i0 + . On the other hand, the pieces containing the integration over the imaginary part of the sunrise remain real since G(t, u) ∈ R if t > u. In order to obtain realvalued polylogarithms it is convenient to perform the change of variables The analytic continuation of the HPLs then is gives The very same steps can be repeated in order to obtain the analytic continuation of the next order, Eq. (9.8). We do not report the results here for conciseness.

Conclusions
The computation of multiloop massive Feynman integrals remains still today an outstanding task due to the appearance of new mathematical structures which cannot be reduced to the by now very well understood multiple polylogarithms. The best known example is that of the two-loop massive sunrise graph. In spite of the recent impressive progress, a formalism which allows to treat not only the sunrise graph, but also, more importantly, more complicated diagrams which, for example, contain it as subgraph, is still missing in the literature. This issue, indeed, becomes of crucial importance for LHC phenomenology, whenever the contribution of massive particles in the loops has to be taken into account.
In this paper we showed that the study of the imaginary part of Feynman graph amplitudes, and the corresponding dispersion relations, can be paired to the differential equations method, providing a very powerful tool for the evaluation of massive Feynman integrals, in particular when the result cannot be written in terms of multiple polylogarithms only. We have considered in detail the case of the kite graph, relevant for the calculation of the two-loop QED corrections to the electron self-energy. The calculation of the kite integral within the differential equations method requires the integration over its full set of subgraphs, which contain both simple integrals which can be expressed in terms of harmonic polylogarithms, and the two-loop massive sunrise. While the former do not constitute any conceptual difficulty and can be treated with standard techniques, the latter require the extension of these techniques. After having established the formalism for the solution of the coupled differential equations satisfied by the two master integrals of the sunrise graph, we showed how to compute their imaginary part and write dispersive relations for the latter. Finally we used these results in order to obtain simple analytical representations for the first two non-zero orders of the kite integral. The final expressions involve polylogarithms up to weight 4 and one-fold integrals over complete elliptic integrals and polylogarithms of weight 2. The numerical evaluation of our result is straightforward, as well as their analytic continuation to all physically relevant values of the momentum squared.
While the problem studied in this paper is relatively simple, the methods presented are very general and can be, in principle, easily extended to consider arbitrarily complicated cases. Moreover, the results derived here, in particular the expressions for the two master integrals of the two-loop massive sunrise, are in a form that is suitable to be re-used once they appear as inhomogeneous terms in the differential equations of more complicated graphs. The application of these techniques to phenomenologically relevant three-and four-point functions is currently under study.
where the four real constants b i satisfy the condition b 1 < b 2 < b 3 < b 4 . We can define three apparently different integrals but in fact they are not all independent. Indeed, consider the contour integral where the contour contains the four points b i . The integrand has two cuts, one cut from b 1 to b 2 , where R 4 (b + i ) = −iR 4 (−b), the other cut from b 3 to b 4 with R 4 (b + i ) = iR 4 (−b). If the contour is the circle at infinity, where 1/R 4 (b) behaves as 1/b 2 , one finds By shrinking the circle to two closed paths containing one of the cuts each, one obtains by comparing the two results for C, one has in general In the case of the equal-mass sunrise the polynomial becomes We define, for n integer and positive, the following three functions such that they are all real-valued as u > 9. Clearly, not all functions are linear independent. Using integration-by-parts identities it is easy to prove that, for each family of functions, only three can be linear independent. We choose for definiteness Moreover one more relation can be written for each family of functions. We find All together these relations imply that only 4 functions are linearly independent. We choose our basis as follows I (0, u), I (2, u), J (0, u), J (2, u).

Appendix B. The analytic continuation of the homogeneous solutions
In the main text we showed how to find the solution of the homogeneous system for the sunrise graph, the matrix G (9,∞) (u) in the region 9 < u < ∞, Eq. (8.10), using the imaginary part of the master integrals as building blocks. In this appendix we show how to build up corresponding real solutions in the remaining three regions, i.e. 1 < u < 9, 0 < u < 1 and −∞ < u < 0.

Appendix D. Relation with the complete elliptic integrals
In this last Appendix we show how to express the solutions entering in the G (a,b) (u) for all four relevant intervals a < u < b, see Eqs. (7.12), (8.10), in terms of the complete elliptic integrals of first and second kind. The latter are defined as (D.1) They are real for 0 < x < 1. From the very definition, one has the particular values where η is small and positive and terms of first order in η are neglected. For η = −ξ − i , with ξ small and positive, > 0 and infinitesimal, the above equation for K(1 − η) gives further K(x), E(x) satisfy the system of first order differential equations given by d dx Considering, more in general, the differential system d dx the pair of functions (K(x), E(x)), obviously provides with a first solution, say F (1) i (x), i = 1, 2 while a simple calculation shows that F (2) 1 (x) F (2) 2 (x) 2 (x) , are also solutions. They cannot be all independent, and in fact one has, for x → x + i , with 0 < x < 1, the relation F (3) i (x) = F (1) i (x) − iF (2) i (x) , i = 1, 2 .
(D. 6) Further, the Wronskian of any two solutions (i, j), 1 (x) is constant (independent of x), as the matrix of the coefficients of the system is traceless. By using the particular values Eq. (D.2) one finds W (1,2) which is the Legendre relation, and