The three-loop QED contributions to the $g-2$ of charged leptons with two internal fermion loops and a class of Kamp\'e de F\'eriet series

The three-loop QED mass-dependent contributions to the $g-2$ of each of the charged leptons with two internal closed fermion loops, sometimes called $A^{(6)}_3\left(\frac{m_1}{m_2}, \frac{m_1}{m_3}\right)$ in the $g-2$ literature, is revisited using the Mellin-Barnes (MB) representation technique. Results for the muon and $\tau$ lepton anomalous magnetic moments $A^{(6)}_{3,\mu}$ and $A^{(6)}_{3,\tau}$, which were known as series expansions in the lepton mass ratios up to the first few terms only, are extended to their exact expressions. The contribution to the anomalous magnetic moment of the electron $A^{(6)}_{3,e}$ is also explicitly given in closed form. In addition to this, we show that the different series representations derived from the MB representation collectively converge for all possible values of the masses. Such unexpected behavior is related to the fact that these series bring into play double hypergeometric series that belong to a class of Kamp\'e de F\'eriet series which we prove to have the same simple convergence and analytic continuation properties as the Appell $F_1$ double hypergeometric series.

m1 m2 , m1 m3 in the g−2 literature, is revisited using the Mellin-Barnes (MB) representation technique. Results for the muon and τ lepton anomalous magnetic moments A (6) 3,µ and A (6) 3,τ , which were known as series expansions in the lepton mass ratios up to the first few terms only, are extended to their exact expressions. The contribution to the anomalous magnetic moment of the electron A (6) 3,e is also explicitly given in closed form. In addition to this, we show that the different series representations derived from the MB representation collectively converge for all possible values of the masses. Such unexpected behavior is related to the fact that these series bring into play double hypergeometric series that belong to a class of Kampé de Fériet series which we prove to have the same simple convergence and analytic continuation properties as the Appell F 1 double hypergeometric series.

Introduction
The Mellin-Barnes (MB) representation method, a well-known computational tool of perturbative quantum field theory, can be used to derive series representations of Feynman diagrams and related quantities in terms of multiple hypergeometric series. In general, once the Feynman diagram or quantity of interest has been represented by a multi-fold MB integral, a standard residue calculation shows that several of such series representations, converging in different parts of the parameter space, can be derived and these series are, as a rule, analytic continuations of one another (see [1] for a systematic exposition of the 2-fold case). However, even at the level of 2-fold MB representations the convergence regions of these analytic continuations do not collectively cover, in general, the whole parameter space of the computed quantity. This implies that one has to find alternative (and sometimes non trivial) analytic continuation methods in order to obtain analytic expressions valid in the particular regions of the parameter space where none of the series derived from the standard residue computation of the MB representation can be used. We call these inaccessible regions the "white regions" in what follows.
A well-known example involving triple series is the two-loop massive sunset Feynman diagram. In [2], two different triple series representations of the latter, derived from its 3-fold MB representation have been given in closed form as combinations of Lauricella F (3) C triple series, and two others can also be obtained, either from the MB representation or by using the invariance of the F parameter space of the sunset diagram (where p is the external momentum and the m i are the masses of the involved particles), but there remains a white region, which includes regions of phenomenological interest, that cannot be reached by any of them. We have shown in [3] how one can analytically continue some of these series to get new series representations of the sunset diagram that can be used to analytically evaluate the latter in several important parts of its white region.
In this paper, we go further on our exploration of the analytic continuation properties of Feynman diagrams and related quantities by revisiting what is possibly the simplest class of QED contributions to the anomalous magnetic moment of each of the charged leptons that can be represented by a 2-fold MB integral. These three-loop QED mass-dependent contributions with two internal closed fermion loops (see Figure 1 for the corresponding Feynman diagram), often denoted A (6) 3 (m 1 /m 2 , m 1 /m 3 ) in the g − 2 literature, can then involve at most double hypergeometric series, and we show that they have an unexpected behavior. Indeed, the analytic continuation properties that these g − 2 contributions satisfy are, surprisingly, the converse of what one faces when one deals with, for instance, the sunset diagram case because the MB representation of A (6) 3 (m 1 /m 2 , m 1 /m 3 ) does not give rise to any white region. This interesting result has encouraged us to probe what is special about the double hypergeometric series involved in our final expressions. Studying the specific form of these series, we observe that they have the same simple convergence and analytic continuation properties as the Appell F 1 double hypergeometric series. Furthermore, we show that the latter and the former both belong to a class of Kampé de Fériet series for which we prove, from their MB representation, the absence of white regions.
As another motivation for studying these particular g − 2 contributions it should be noted, and as emphasized in [4], that in contrast to all other three loop QED contributions to the muon anomalous magnetic moment, A 3,µ (m µ /m τ , m µ /m e ) is the only one whose exact × l 1 l 2 l 3 Figure 1: The 3-loop QED Feynman diagram corresponding to A analytic form has not been derived so far. Results were first presented in [5] in terms of the first few terms of a series expansion in powers and logarithms of the mass ratios, using largemomentum, heavy mass and eikonal expansions techniques. These results have then been checked and extended in [6] using the MB representation method. In the present paper, we have derived them in their entirety and present their exact expressions, in terms of generalised hypergeometric and Kampé de Fériet double hypergeometric series. In the case of the electron, we have not been able to find any analytic result for these contributions in the g −2 literature, although some numerical evaluations of these have been given, for instance in [7]. We will show in the following that the exact expression of A 3,e (m e /m τ , m e /m µ ) has a simple and compact form. The τ lepton case is more intricate and has been considered a long time ago in [8]. The latter reference gives, to our knowledge, the only available non-numerical result for A (6) 3,τ (m τ /m µ , m τ /m e ). The result of [8] corresponds to the leading term in the double series expansion of the exact expression which we will present in the following. The numerical evaluation of this leading term, presented in [9], does not agree with the numerical evaluation of A (6) 3,τ (m τ /m µ , m τ /m e ) given in [7]. We show here that this mismatch can be solved once one adds some sub-leading terms to the expression of [8,9].
In view of all the considerations spelt out in the foregoing, we now give the outline of this paper. In Section 2, a short review of the QED contributions to the anomalous magnetic moment of charged leptons is given. In Section 3 we present the MB representation for the three-loop contribution to g − 2 coming from the Feynman diagrams of Figure 1, and then calculate it for the cases with external electron, muon and τ lepton. Detailed expressions for each of these are listed in the Appendix. In Section 4, we present the checks of our formulas and give a brief numerical analysis, and in Section 5 we discuss the analytic continuation properties of the class of Kampé de Fériet series mentioned above. We conclude with Section 6, where a short discussion of the results and future work are presented.

Short QED literature review
The anomalous magnetic moment of the charged leptons is defined as a l ≡ (g l − 2)/2, where g l is the Landé factor and l = e, µ, τ . In the Standard Model, contributions to a l arise from electroweak and strong processes. The anomalous magnetic moment of charged leptons has a distinguished place in elementary particle physics. Historically, the electron anomalous magnetic moment has been among the most important tests of quantum electrodynamics (QED). During the last decade, a persistent discrepancy between the Standard Model theoretical pre-dictions and experimental results in the case of the muon, now reaching a 3.5 σ level, has spurred a new experiment at Fermilab, from which a 5 σ deviation could well be obtained in the near future [10]. Due to these reasons, over the decades a huge amount of theoretical effort has been devoted to computing this quantity (see [4] for a recent and comprehensive review in the muon case), with a great deal of activity being directed at computing hadronic contributions to it. We recall that, due to the larger mass of the muon, it is usually accepted that the muon anomalous magnetic moment is more sensitive to new physics than the electron. The very short life time of the τ lepton prevents the experimental measurement of its anomalous magnetic moment, which explains why the theoretical study of the latter is less well developed (see [11] for a review of the τ lepton case). The best experimental limits are −0.052 < a τ < 0.013 (95%CL) [12].
The QED contributions to the anomalous magnetic moment of the charged leptons can be expressed perturbatively as is the sum of the n th loop contributions, and A 3,l = 0. The A 1,l are mass independent and thus equivalent for all three lepton flavours. Up to O(α 3 ), these contributions are known in closed analytic form, while A and higher-loop contributions are only known numerically. That A (2) 1 = 1/2 has long been known [13], as has the value of A 1 [14,15,16]. A numerical value of A (6) 1 is given in [17], and after work spanning several decades in calculating the various diagrams that constitute A (6) 1 (e.g. [18,19,20,21,22,23,24]), the calculation was finalised and the exact close form analytic result was presented in [25]. Purely numerical calculations for A (8) 1 can be found in [26,27,28,29,30,31,32], and a result to 1100 digit precision, accompanied by a semi-analytic fit to the result, is presented in [33]. A numerical value for A (10) 1 is given in [31,32]. The mass dependent terms A (4) 2,l have been calculated as a series expansion in the small mass ratio [34,35,36,5], as an exact result for small mass ratios [37], and finally as a closed form analytic result for all values of the mass ratio [38]. As the diagrams under consideration in this work can contribute at the three-loop level when two of the three involved leptons have the same flavor, we discuss A (6) 2,l in greater detail below. The fourth and higher loop results were known primarily numerically, e.g. [39,31,32], although some analytic results are also available, e.g. [40,41,42,43,44].
One of the earliest calculations of A 2,µ was done in [45], where an expansion including some of the leading log and analytic terms was presented. This expansion was extended in [46], and a result for Figure 1, with l 1 = µ and l 2 = l 3 = e was given with log contributions up to order (m µ /m e ) 0 included. In [47,48,49,18] calculations and results for some diagrams making up A (6) 2,µ and A (6) 2,e are presented. In [50] analytic results for vacuum polarization contributions to A (6) 2,µ and A (6) 2,e , up to order (m µ /m e ) 0 , are given. Expressions for all the individual graphs (such as for Figure 1) are however not presented. [51] completes the calculation of all the log(m µ /m e ) terms by computing light-by-light scattering diagrams, numerical estimates for which are to be found in [52,53]. The electron light-by-light scattering graph contributions to A (6) 2,µ are calculated in [54], and the expansion presented there is extended to higher order in [55]. [35] continues the expansion given in [50] to order (m µ /m e ) 0 . In [56] close form analytic results for the vacuum polarization diagrams contributing to A (6) 2,µ are given, as are expansions in mass ratios (up to certain order) that can be used to calculate the equivalent contributions to A (6) 2,e and A (6) 2,τ . The exact expression given in [56] is also expanded in the small mass ratio in [57], but to a higher order than in the former, and which is then used to calculate and compare with numerically evaluated values of A (6) 2,e . The three loop contribution that we are primarily concerned with in this work is the lowest order non-zero term consisting of three masses (or two mass ratios), A 3,l . In [35], an integral representation is given for A (6) 3,l and numerically evaluated for l = µ (we have noted that Eqs. (32) and (33) of [35] do not match numerically because of a 2 overall factor that is missing in the r.h.s. of Eq.(32)). Another numerical evaluation of A (6) 3,µ is given in [56]. For A (6) 3,µ , [5] gives an expansion up to the first few terms, based on asymptotic and eikonal methods, while in [6] these results are slightly extended, using the Mellin-Barnes technique. In [39], results for QED contributions to muon g − 2 up to the tenth order are given, which includes numerical values for A (6) 3,µ and A (6) 2,µ . The integral of [35] is evaluated for the electron case, and a numerical value for A (6) 3,e is given in [57]. Numerical results for the same, as well as for A (6) 2,e , but based on an older set of mass value inputs, is given in [41]. And [31] gives numerical results for QED contributions to electron g − 2 up to the tenth order, including values for A (6) 3,e and A (6) 2,e . [9] reviews and updates all contributions to the muon and τ lepton g − 2 up to its date of publication, and based on the results of [8], gives an expression for A (6) 3,τ as an expansion up to its leading term, which corresponds to our R {3,1} result, see Section 3.3, Table 5. Its numerical value is 2.75316 [9] and can be compared to the numerical value for A (6) 3,τ given in [58]: 1.679. The numerical result for A (6) 3,τ calculated in [57] (see also [11]) is 3.34797, which therefore disagrees with those of [9,58] but agrees with our own calculated value as shown later on in this work. Since the result of [58] is precisely half of the one of [57], we suppose that, as mentioned above, the factor of 2 missing in the r.h.s. of Eq.(32) of [35] is at the origin of this discrepancy, since the authors of [35] are also those of [58]. Concerning the discrepancy with [8,9], as already said in the introduction, adding a few sub-leading terms allows to obtain agreement with [57].
For a recent review of the theory and experimental status of the g − 2, see [10] ( [59] is also useful). For a comprehensive review of QED contributions to all the leptons' g − 2, see [60,57]. A review of the muon g −2 is given in [61,4,62,63,38], and a review of contributions to τ lepton g − 2 is presented in [11].

Three-loop QED contributions with two internal loops
We will now give the exact expressions of the contributions A (6) 3,l to the anomalous magnetic moments of the electron, muon and τ lepton coming from the Feynman diagram of Figure 1 (and the symmetric diagram obtained by an exchange of the internal loops).
The MB representation of these contributions may be found in [6]. Defining r 1 .
where γ . Figure 2). From the rules described in [1,64] it is clear that we are in a so-called degenerate case (∆ = 0), where several convergent series representations of the integral coexist, being, as a rule, analytic continuations of one another.
Since the MB integral is fully symmetric under the exchange of r 1 and r 2 (or s and t) one can avoid about half of the calculations that would be necessary to perform in order to derive all the possible convergent series representations in the case of a non symmetrical integral. This symmetry, which comes from the symmetry of the Feynman diagram under the exchange of l 2 and l 3 , is also reflected in the singular structure of the integrand (see Figure 2) and in the picture showing the convergence regions of the series representations in the first quadrant of the (r 1 , r 2 )-plane (see Figure 4). It is easy to find the different sets of residues (cones) associated to each convergent series (see [1] for details on the general procedure). There are six such cones, plotted in Figure 3 and, as just explained, the series representations associated to only three of them have to be computed (the blue cones), the others (coming from the red cones) being obviously derived from the latter by exchanging r 1 and r 2 in the final results.
We show in Figure 4 the convergence regions of the series representations deduced from each of the cones and, in Table 1, we show to which of the cones the possible physical situations are associated. The mass independent case can be computed from any of the cones, as can be seen in Figure 4 .

The electron case (Cone 1)
Let us begin with the simplest case, namely A 3,e : a muon and a tau in the internal loops and an electron on the external legs of the Feynman diagram shown in Figure 1, and its symmetric counterpart. We will see that the exact analytic expression of these contributions to the anomalous magnetic moment of the electron is more compact than the corresponding expressions for the muon and tau cases. As mentioned in the introduction, we have not been able to find any analytic result for this contribution in the g − 2 literature.
A convergence analysis to be presented below shows that the series representation corresponding to the electron case comes from Cone 1 (or Cone 4 for the symmetric diagram) in m τ m µ m µ 3,6 Table 1: The different combinations of charged leptons and corresponding cones. l 1 is the external lepton (see Figure 1). associated poles in Figure 3 are : • Single series contributions: • Double series contributions:

Exact result
The series representation of these contributions to the anomalous magnetic moment of the electron is expressed in closed form in terms of the generalized hypergeometric series 4 F 3 and the Kampé de Fériet double hypergeometric series F 2:3;3 2:2;2 and F 2:3;2 2:2;1 , and the final expression reads The residues R {1,i} are given in Table 2 as well as their correspondence to the singular points listed above. As the expressions for some of the residues is lengthy, we have introduced some functions, h {1,j} , in the table, whose explicit forms are given in the Appendix, see Eqs.
Due to the presence of non-simple poles in the singularity structure of Eq.(3), many of its residues involve polygamma functions, which come from derivatives of the gamma function. As these polygamma terms arise from an application of the residue theorem, they are derivatives of gamma functions appearing in the h {i,j} of the concerned singularity. Therefore, it is always possible to express the residues solely in terms of gamma functions, and to express the polygamma factors as derivatives of those gamma functions. As an example, the explicit Singularity Label Residue The above was obtained by applying the following Cauchy's theorem operator and thereafter setting s = t = 0 (see [1] for details on the general computational procedure). R {1,2} can therefore be expressed more compactly as where and where we have replaced the variables s and t by the parameters α and β, and expressed the sum over m in terms of the generalised hypergeometric function, 5 F 4 . The advantage of Singularity Label Residue this notation is that it is concise, and that by expressing the single series as p F q , or the double series as Kampé de Fériet series one may perform analytic continuations on these results easily if needed. It is in some cases possible to further simplify the results and express them in terms of a single parameter. For example, R {1,2} may be expressed as where In the rest of this paper, we have chosen to express the results in the most compact notation possible. However, for illustrative purposes, we have given both forms (i.e simplified and non-simplified) of the results in the electron case (see Table 2 and Table 3).
The h {1,i} and the H {1,i} are given in the Appendix.

Convergence region: external electron
The region of convergence of the corresponding series representation is straightforward to derive from the convergence properties of generalized hypergeometric and Kampé de Fériet series [65]. For the generalized hypergeometric series they read (assuming that none of the parameters is zero or a negative integer): then the p F q series, with p = q + 1 is For the Kampé de Fériet double hypergeometric series (with the following notation in the l.h.s of Eq. (14): (a p ) . = (a 1 , ..., a p )) we have (i) Convergence for |x| < ∞ and |y| < ∞ if p + q < l + m + 1 and p + k < l + n + 1.
(ii) If p + q = l + m + 1 and p + k < l + n + 1, the convergence is for |x| From the results above and Table 2 it is easy to find that the convergence region of the r.h.s. of Eq.(4) is simply This region is plotted in Figure 4

The muon case (Cone 2)
The muon case corresponds to Cone 2 in Figure 3. Therefore, the different sets of singularities to consider in the (Re(s), Re(t))-plane are the following, where as before m and n are any non-negative integers: Three of them have already been considered during the calculation of Cone 1, in the electron case. One should however keep in mind that in Cone 2, their corresponding transformation law [1] will not be the same as in Cone 1 so that these three sets of residues will not give the same analytic expression in Cone 1 and Cone 2.

Singularity Label
Residue

Full result
The series representation extracted from Eq.(3) by summing the residues of Cone 2 is: As in the electron case, the correspondence of the singularity points and their residues R {2,i} are presented in Table 4, the explicit forms of some of the residues being relegated to the appendix for lack of space in the main body of the paper (see Eqs.(B-1)-(B-8)).

Convergence region: external muon
Using the results presented in Section 3.1.2, it is easy to conclude that the convergence region of the r.h.s. of Eq. (16) is See Figure 4 for a plot of this region (labelled "Cone 2"). As in the electron case, it is possible to include the boundaries in this convergence region. In the case where the three different leptons are involved, the only phenomenological situation which satisfies Eq.

The τ lepton case (Cone 3)
The τ lepton case falls in the convergence region associated to Cone 3 (see Figure 3), and it is the hardest from the computational point of view for two reasons: there are a lot of different types of singularities in this cone, and here one also has to take care of the cancellation, or reduction of multiplicity, of different sets of singularities due to gamma functions in the denominator of the MB integrand (a cancellation also happened in the muon case but only for one set of singularities). With m and n any positive integer, one may as usual exhibit all sets of singularities contributing to the cone. There are 25 different sets to consider:

Full result
The series representation extracted from Eq.(3) by summing the residues of Cone 3 is: As in the preceding electron and muon cases, the correspondence of the singularity sets and their residues R {3,i} are presented in Table 5 and the explicit forms of some of the residues is relegated to the appendix (see Eqs.(C-1)-(C-15)).

Convergence region: external tau lepton
Once more, using the results presented in Section 3.1.2, one concludes that the convergence region where the series representation associated to Cone 3 is valid is: See Figure 4 the region labelled "Cone 3". As in the two previous cases, it is in fact possible to include the boundaries in this convergence region. Here we see that with r 1 = m 2

Other cones
As already mentioned, due to the symmetry of the MB integral, it is possible to obtain the results of the other three cones by a simple interchange of r 1 and r 2 in the results that we have already obtained. One then sees in Figure 4 that the whole first quadrant of the (r 1 , r 2 )-plane may be reached. Values at the boundaries of the different cones may be evaluated by using expressions of either cone. This is discussed further in Section 4 and 5. This behaviour is therefore completely different from what can be found in the examples considered in, for instance, [1,2,66,67] where there were always white regions in the parameter space which could not be reached.

Singularity Label
Residue We performed several internal and external checks to ensure the validity of our expressions. The first consistency check was to compare values obtained by a high precision numerical integration of Eq.(3) and its Feynman parametrization against our analytic results. For each cone, we have tested our expressions by ensuring that there is agreement between the integral and our full analytic result, both computed numerically, to at least six orders of magnitude beyond the order of magnitude of the smallest contributing set of residues of that cone.
The second consistency check was to compare our results for Cone 2 and Cone 3 against another set of analytic expressions derived by analytically continuing the residues of Cone 1, i.e. R {1,i} for i = 1, ..., 4. This results in a set of series that are numerically equivalent to R {2,i} and R {3,i} , but which are different in their analytic form. These were numerically evaluated, and shown to agree with the A 3,µ and A 3,τ computed directly from the R {2,i} and R {3,i} , respectively, to at least 19 decimal places. Let us briefly describe this approach on the example of the first electron residue R {1,1} (see Table 2), which involves a Kampé de Fériet series. The latter has the following Mellin-Barnes type integral representations [69] and . Solving this MB representation in appropriate cones yields the desired analytic continuations for that particular series.
The convergence regions of the series expansions obtainable from a direct residue calculation of Eq. (22) following the method of [1] are shown in Figure 5 and labelled by their associated cone. (Note that we use roman numerals to label and distinguish the cones of Eq. (22) from those of Eq.(3)). The sum of residues of Cone i obviously reproduce the Kampé de Fériet series of R {1,1} . For the muon case, where r 1 = m 2 µ /m 2 τ ∼ 10 −3 and r 2 = m 2 µ /m 2 e ∼ 10 4 , the residues of Cone iv have to be calculated. For the tau case where r 1 = m 2 τ /m 2 µ ∼ 10 2 and r 2 = m 2 τ /m 2 e ∼ 10 7 , those of Cone v need to be calculated.
By performing a similar calculation on R {1,2} , R {1,3} and R {1,4} , and summing the results of the appropriate cones, the results of Cone 1 (external electron legs) of Eq.(3) may be analytically continued to obtain the results of Cone 2 (external muon legs) and Cone 3 (external tau legs) of the same integral; see [70] for details and results of the complete calculation. This results in series representations that are numerically equivalent, but different in form, to those obtained from a direct calculation of Eq.(3). This process did not yield expressions that were simpler than those obtained by a direct evaluation. However, results obtained by this analytic continuation approach could well confer advantages, such as simplicity of form, in other cases than the g − 2. We also checked our results externally by numerically comparing them with results from the literature.
For the muon, our expressions evaluated with the CODATA 2010 [71] mass ratios yield the same value of A (6) 3,µ as given in Eq.(9) of [39]. Similarly, we get agreement with Eq.(22) of [56] evaluated using the values of the 1992 PDG [72], and with Eq.(16) of [38] evaluated using CODATA 2002 [73] inputs, the latter of which is based on the first few terms of the asymptotic and eikonal expansion derived expressions of [5]. Our expression differs numerically from the value given in Eq.(33) of [35] at the sixth decimal place. The same integral is evaluated in [57,11] with CODATA 2006 [74] mass ratios to yield A (6) 3,µ , values with which our expressions are in agreement. Our results for A (6) 3,τ , however, yield a numerical value that differs from the older literature values [9,58] but agrees with [57], as already mentioned in the introduction and in Section 2. For the electron, our A (6) 3,e expression agrees numerically with the value given in [57], based on the integral of [35] and 2004 CODATA values, as well as with the result given in Eq.(12) of [41] that uses the PDG 2012 values as inputs. And finally, within the latter's uncertainty range, we obtain the same result as [31] using the CODATA 2010 [71] mass ratios.
Note that the diagram of Figure 1 appears in the expansion of a l not only as A 3,l , but also as one of the several diagrams that constitute A  Table 1 gives, among others, the various mass configurations possible for Figure 1 with two masses, and the cones corresponding to their evaluation. As all these cases lie on the boundary of the convergence regions of two cones (in passing, the mass independent case can be computed with any of the 6 different series representations) the convergence of the series can be slow. However, by comparing the literature values of these diagrams with our expressions, we obtain another check of the latter. Note that for the cases where there are identical leptons in the loops, an overall factor of 1/2 is required in front of Eq. 3.
By setting all masses equal, we obtain the scaleless diagram A 1,l , which corresponds to the 'triple point' (1.0, 1.0) in Figure 4. Despite the slow convergence of our series at this point, we find that as the number of terms in the sums is increased, the results of all cones numerically tend to the value given by Eq.(8) or (9) of [56] or more explicitly by Eq.(4.43) of [43]: 8 3 The same equations of the paper [56] give the closed form result for the cases where both leptons in the internal loops are identical, but different from the one on the external legs, and where only one of the leptons in the internal loops is the same as the one on the external legs. And Eqs. (12), (13), (16) and (17) of [56] give the expansions in the mass ratio of Eqs. (8) and (9). For the case where the muon is on the external legs, with two internal electron loops, we get from our expression a fast converging numerical agreement with the value calculated using Eq.(8) of [56]. With two internal tau loops the agreement is also very good, but we use Eq.(16) of [56] to compare because imaginary contributions appear when using Eq. (8). With only one internal electron or tau loop, we get values that are converging in a slower way to those obtained using Eq.(9) of the aforementioned paper as we increase the number of terms in the sum. The same is true in case of tau external legs, where we get very good and fast converging numerical agreement with Eq.(8) or (9) of [56] for loops with two electrons or two muons, and slower convergence towards the literature values otherwise. For the case with electron masses in the external legs, we do not use the closed form expressions of Eq.(8) and (9) of [56], but rather the expansion in the mass ratios given in Eqs.(16)-(17) of the same paper. As with the other cases, we get very good agreement when the two fermionic loops carry the same mass, and slower convergence otherwise. All these "two masses" results are also in agreement with the numerical integration of the Mellin-Barnes representations given in Eqs.(4.23) and (4.25) of [43]. cover the entire parameter space (except possibly when the absolute value of the parameter is equal to unity).
In our present case of study, where N = 2, it is then clear that the lack of white regions is related to the particular properties of the double series that appear in the series representations derived in the previous sections. A quick look in the Appendix shows that there are two different types of double series in all the mathematical expressions derived from the different cones, namely the Kampé de Fériet series These Kampé de Fériet double series also appear in our formulas with different arguments, but one can perform the analysis of this section to these cases also and the conclusions will be the same. Note that Eq.(24) can in fact trivially be transformed to the form of Eq.(23) by including a (c 3 )n (c 3 )n factor in the sum over n of the definition of Eq.(24) given in Eq. (14). Therefore, let us focus on Eq.(23).
From Section 3.1.2, one can deduce that the convergence region of F 2:3;3 2:2;2 is the simple region |x| < 1 ∧ |y| < 1. And in Section 4 we have seen that from the MB representation of F 2:3;3 2:2;2 (given in Eq. (22) for specific values of its coefficients, which have no consequence on the convergence properties), one can show that the convergence regions associated to each of the possible series representations derived from the MB integral are those of Figure 5, where there is no white region.
In fact, the convergence and analytic continuation properties of the Kampé de Fériet F 2:3;3 2:2;2 double hypergeometric series are exactly the same as those of the Appell F 1 double hypergeometric series. It is easy to show this from the MB representation of the latter. The similarity between the Appell F 1 and the Kampé de Fériet F 2:3;3 2:2;2 series should not come as a complete surprise because Kampé de Fériet series are generalisations of Appell series in the same way as generalised p F q hypergeometric series are generalisations of the Gauss 2 F 1 hypergeometric series. Obviously, in the latter case there is only one function to generalise, whereas in the Appell case, there are four different ones which do not satisfy the same convergence and analytic continuation properties. In particular, F 1 is the only Appell series whose MB integral representation gives birth to some series representations, analytic continuations of one another, whose convergence regions are able to collectively cover the whole (|x|, |y|)-space (except possibly on some boundary lines at which the convergence depends on the values of the coefficients). Indeed, for the other Appell and Eqs. (23) and (24), it is natural to conjecture that the more general Kampé de Fériet series of the type x, y where we recall that (a 1+i ) . = (a 1 , ..., a 1+i ). To simplify our presentation, we suppose that the a j , b j and c j are non-zero positive numbers. This allows us to deal with straight contours and γ . = (Re(s), Re(t)) therefore belongs to the non empty fundamental polygon defined by the positivity constraint of each of the real parts of the gamma functions in the integrand. We also suppose that the singular structure of the integrand has no poles of multiplicity greater than 1, and that the values of the α j , β j and γ j cannot make the gamma functions of the denominator interfere with those of the numerator.
Following the method of [1] it can be shown that this MB integral has five different cones and therefore five different double series representations. We list the sets of singularities for each cone in Table 6 (m and n can be any positive integer or zero).
In the g − 2 case, the MB integral of Eq.(3) is symmetric under the exchange of r 1 and r 2 , and therefore of the total number of series, three of them could be trivially obtained from the three others using that symmetry. Here, one can still use this symmetry, although since k = l and b j = c j generally in Eq. (27) one has to exchange, in addition to r 1 and r 2 , the values of the coefficients b j and c j as well as k and l. Since Cone 1 obviously gives the Kampé de Fériet series in Eq. (26), it is sufficient to compute the series representations associated with Cone 2 and Cone 3, the two others being deduced from them.
Let us begin with Cone 2 which has only two different types of contributions (see Table  6). In this case, one obtains where the * superscript means that, for instance, the r = j case is not considered in the product 1+i, * r=1 Γ(a r − a j ), and that ((1 − a + a j ) 1+i ) * = (1 − a 1 + a j , 1 − a 2 + a j , ..., 1 − a j−1 + a j , 1 − a j+1 + a j , ..., 1 − a 1+i + a j ).
It is easy to see, in the analytic continuation formula Eq. (28), that the first series, being a Kampé de Fériet double hypergeometric series, converges in the region y x < 1 ∧ 1 x < 1. The second series in Eq. (28) is not a Kampé de Fériet series, but by the cancellation of parameters method [65], one finds that it converges as a Horn G 2 double hypergeometric series, i.e in the region 1 x < 1 ∧ |y| < 1, so that the convergence region of the analytic continuation of F 1+i:1+k;1+l 1+i:k;l associated to Cone 2 and given in Eq.(28) converges for 1 x < 1 ∧ |y| < 1. This region corresponds to the green region plotted in Figure 5.
One can now proceed to the presentation of the results associated with Cone 3. In this case, one gets A convergence study similar to the one performed above shows that one finally obtains, for the series representation presented in Eq. (29), the same convergence region as the one plotted in orange in Figure 5, i.e the region y x < 1 ∧ 1 x < 1 ∧ 1 y < 1. It is straightforward to show that the series associated with Cone 4 and Cone 5 will converge, respectively, in the purple and in pink regions of Figure 5, and one can check that our conjecture that there are no white regions for the class of Kampé de Fériet series considered in Eq.(26) is correct. Therefore, we see that when dealing with this class of series, one does not need to perform non trivial analytic continuations in white regions, such as is needed for those of [3] (see also [75]).
We conclude this section by noting that a similar analysis can be performed in order to extract other classes of Kampé de Fériet series that have common convergence and analytic continuation properties to the three other Appell series (and therefore also the same white regions). We do not perform such an analysis here but instead give two examples of Feynman diagrams taken from some of our recent work where there are similarities between the Kampé de Fériet series found in their analytic expressions and the Appell F 3 and F 4 series. The two-loop box diagonal calculation presented in [76] involves a Kampé de Fériet F 2:2;2 3:0;0 double series which has the same white region as the Appell F 3 series. Let us recall that The example related to the Appell F 4 series comes from the chiral perturbation theory sunsets studied in [66]. In theH π KKη pion andH η KKπ eta sunsets analytic expressions, one can find the Kampé de Fériet F 3:1;1 1:2;2 double hypergeometric series, which has the same white region as These two examples related to F 3 and F 4 need analytic continuation procedures alternative to the traditional MB representation in order to be analytically computed using convergent series representations in their white regions.
The discussion of the present section draws attention to the fact that the "Appell F 1 " class of Kampé de Fériet double hypergeometric series presented in Eq.(26) has a particularly simple analytic continuation behaviour and it suggests that these nice properties are also very probably satisfied by extensions of Kampé de Fériet series. Such a higher order class of Kampé de Fériet multiple series would then include the F (n) D Lauricella series as the simplest series of its set.

Conclusion and Discussion
The aim of this work was two-fold: to present complete analytic results for the three loop contributions to leptonic g-2 with two internal loops, and to use these calculations to further our understanding of the convergence and analytic continuation properties of multiple hypergeometric series that can appear in quantum field theory calculations.
In the first part of this work, we have calculated and presented complete analytic results for the three-loop QED contributions to the g − 2 of all charged leptons with two internal loops, i.e. for all three possibilities of external legs in the Feynman diagram of Figure 1, denoted in the literature by A (6) 3,l for l = e, µ, τ . In the muon case, this was the last missing piece in the puzzle of the exact results at three loop level [4]. Furthermore, to our knowledge, in the electron and τ lepton cases, analytic results were unknown but for the leading term in a double expansion in the mass ratios for the τ lepton case given in [8,9]. Therefore, this work presents the first complete and exact analytic result for A (6) 3,l for all l = e, µ, τ . We have performed several checks of the expressions given in this paper. These included numerical comparisons with values from the literature, as well as consistency checks. One such consistency check involved the calculation of the diagram of Fig. 1 with two or three of the leptons being identical, i.e. the contributions of Fig. 1 to A (6) 1,l and A (6) 2,l . The use of only two distinct masses, or one single mass, in expressions consisting of two mass ratios lands us on the boundaries of the convergence regions, which may be reached by two or more distinct series. Finding numerical agreement for expressions that were different in form was one consistency check (note that for these cases we have also found agreement with numerical results derived from the analytic expressions given in [56] and from the numerical integration of the Mellin-Barnes representations of [43]). The same principle was the basis for our other self-consistency check. In this check, we have calculated a second set of expressions by analytic continuation of the results obtained in the electron case (see [70] for a complete listing of these expressions), and then compared them to Eq. (16) and Eq. (18), which are numerically equivalent but have a different form.
Our tool for performing the calculation was the Mellin-Barnes representation, which produced results in the form of a linear combination of isolated terms and infinite single and double series consisting of products of gamma and polygamma functions, that can then be expressed in terms of generalized hypergeometric and Kampé de Fériet series and their derivatives. These series in turn become the objects of study of the second part of this paper.
The second part of our work involved studying the analytic continuation of a class of Kampé de Fériet series. Indeed, the generalized hypergeometric and Kampé de Fériet series converge in a range of values of their variables. By analytic continuation, one is able to extend this range. But for the double hypergeometric series, when this analytic continuation is performed by means of the series derived from standard residues computations of MB representation, there usually exists a range of values (the white region) for which it is still not possible to derive some converging series. However, for the series appearing in the g − 2 calculation of this paper, we find no white regions.
Inspired by this unusual analytic continuation property, in this paper we prove that for a class of Kampé de Fériet series, namely those of type F 1+i:1+k;1+l 1+i:k;l , no white region can appear, and that it is possible to analytically continue these series using their MB representation to obtain convergent series for all values of their variables (except possibly, depending on the values of the parameters, on the boundaries of the convergence regions). The convergence and analytic continuation properties of the F 1+i:1+k;1+l 1+i:k;l series parallel those of the Appell F 1 series. We also give examples of some Feynman diagrams that indicate that other classes of the Kampé de Fériet series, whose convergence and analytic continuation properties mimic those of the other Appell series, may be found in physical situations.
In this work, therefore, we present complete analytic results at the three-loop level for the A 3,l contribution to the important physical quantity g − 2, and in the process of the calculation extend our understanding of the convergence and analytic continuation properties of the Kampé de Fériet series. Further investigating these properties for other classes of Kampé de Fériet series and multiple series of higher order, of which relatively little is known, is an important direction for future research.