The Mellin-Barnes Approach to Hadronic Vacuum Polarization and $g_{\mu}-2$

It is shown that with a precise determination of a few derivatives of the hadronic vacuum polarization (HVP) self-energy function $\Pi(Q^2)$ at $Q^2=0$, from lattice QCD (LQCD) or from a dedicated low-energy experiment, one can obtain an evaluation of the lowest order HVP contribution to the anomalous magnetic moment of the muon $a_{\mu}^{\rm HVP}$ with an accuracy comparable to the one reached using the $e^+ e^-$ annihilation cross-section into hadrons. The technique of Mellin-Barnes approximants (MBa) that we propose is illustrated in detail with the example of the two loop vacuum polarization function in QED. We then apply it to the first few moments of the hadronic spectral function obtained from experiment and show that the resulting MBa evaluations of $a_{\mu}^{\rm HVP}$ converge very quickly to the full experimental determination.

(1. 5) We observe that the integrand in Eq. (1.4) can be rearranged in a way: which explicitly displays the dispersion relation between the hadronic spectral function and the renormalized hadronic photon self-energy in the euclidean: and therefore [21,22] a HVP (1.8) Trading the Feynman parameter x-integration by a Q 2 -integration results in a slightly more complicated expression which is the one proposed for LQCD evaluations [23]. Because of the parametric x-dependence in Eq. (1.8), or the Q 2 -weight function in the integrand of Eq. (1.9), the a HVP µ integral is dominated by the low-Q 2 behaviour of the hadronic self-energy function Π(Q 2 ). The natural question which then arises is: What is the best way to help LQCD (see e.g. refs. [10]- [19]), or dedicated experiments [20], to evaluate this integral when only limited information about Π(Q 2 ) at low Q 2 values is available? The answer that we propose follows the way initiated in ref. [2]. It is based on Mellin-Barnes techniques which we shall describe below and which we shall illustrate with several examples. As we shall see, this is a very powerful method compared to other approaches discussed in the literature (see e.g. refs. [24,25,27] and references therein).
The paper has been organized as follows. The next section is an introduction to the QCD properties of the Mellin transform of the HVP spectral function. Section III is dedicated to a few ingredients, which are required to understand and justify the method that we propose. The subsection III.3 is particularly technical since it justifies mathematically the underlying approach and the restriction to the subclass of Marichev-like Mellin approximants given in Eq. (3.16). For those who are just interested in the applications, it can be escaped in a first reading. Section IV illustrates the application of Mellin-Barnes approximants (MBa) to vacuum polarization in QED at the two loop level. Section V tests the advocated technique of MBa with the experimental values of the HVP moments provided to us by the authors of ref. [9]. These moments, with their errors, are obtained from the same spectral function which results in the second number quoted in Eq. (1.3). We show how the successive MBa approach the experimental determination of a HVP µ . The conclusions with an outlook on future work are given in Section VI. A few technical details have been included in an Appendix.

II The Mellin Transform of the Hadronic Spectral Function.
In QCD the hadronic spectral function is positive and goes asymptotically to a constant (q i denotes the charge, in electric charge units, of an active quark with flavour i ) : π ImΠ(t) , n = 0, 1, 2 · · · , (2.2) where throughout the paper t 0 denotes the threshold value of the hadronic spectral function: can be experimentally determined; and the dispersion relation in Eq. (1.7) relates them to successive derivatives of the hadronic self-energy function Π(Q 2 ) at the origin: , n = 0, 1, 2, · · · , (2.4) which are accessible to LQCD evaluations. In fact, as pointed out a long time ago [26], the first moment for n = 0 provides a rigorous upper bound to the muon anomaly: . (2.5) Quite generally, the moments in Eq. (2.2) obey constraints which follow from the positivity of the spectral function and may provide useful tests to LQCD determinations. We discuss these constraints in the Appendix. The moment integrals in Eq. (2.2) can be generalized to a function, which is precisely the Mellin transform of the hadronic spectral function 1 π ImΠ(t) defined as follows [1]: As a result, M(s) can have neither poles nor zeros in the negative Re(s) axis and has a perfectly smooth (increasing) shape in this region. This smoothness property of M(s), which is at the basis of the approximation method that we shall propose, is to be contrasted with the shape of the spectral function 1 π ImΠ(t) itself which, as we know from experiments, has a rather complicated structure. In QCD, the Mellin transform M(s) is singular at s = 1 with a residue which is fixed by the pQCD asymptotic behaviour of the spectral function in Eq. (2.1). The contribution from the u, d, s, c, b and t quarks gives M(s) ∼ s→ 1 α π 4 9 + 1 9 + 1 9 + 4 9 + 1 9 + 4 9 N c 1 3 The spectral function moments are, therefore, the particular values of the M(s) function at s = 0 , −1 , −2 , −N with integer N . As discussed in refs. [1,2] there exists a representation of Π(Q 2 ), and hence of the anomaly a HVP µ , in terms of the Mellin transform M(s). This follows from inserting the Mellin-Barnes identity 1 in the dispersion relation in Eq. (1.7), which results in the representation 2.10) and the corresponding integral representation for the Adler function The integral over the x-parameter can now be made analytically, leading to the expression [1] a HVP where F(s) is a product of three Gamma functions: 15) and the hadronic dynamics is thus entirely factorized in the Mellin transform M(s).
The weight function F(s) in Eq. (2.14) is universal and has a shape which, for s within the fundamental strip [28]: c s ≡ Re(s) ∈]0, 1[ and the choice s = 1 2 − iτ , is shown in Fig. (1) as a function of τ . Notice that the real part of this function (the red curve) is symmetric under τ → −τ while its imaginary part is antisymmetric. Both the real and imaginary parts fall very fast as τ increases. With the change of variable the integral in Eq. (2.14) becomes then a Fourier transform: 1 For the benefit of the reader who may be unfamiliar with Mellin-Barnes integrals we give a proof of this identity in the Appendix. Shape of the function F 1 2 − iτ in Eq. (2.14) versus τ . The red curve is the real part of the function, the blue dashed curve its imaginary part.
Because of the shape of the F 1 2 − iτ function and the growth restrictions on M 1 2 − iτ for large τ , which are fixed by the fact that Π(Q 2 ) obeys a dispersion relation in QCD, this Fourier integral is fully dominated by the behaviour of the integrand in a very restricted τ -interval, −T ≤ τ ≤ +T with T of order one.
III Some Technical Ingredients.
We shall next recall a few technical ingredients which in the literature go under the name of: Ramanujan Master Theorem, Marichev class of Mellin transforms, Generalized Hypergeometric Functions and Meijer's G-Functions. They are necessary to implement and justify the MBa framework that we propose.

III.1
The so called Ramanujan's Master Theorem.
Consider a function F (x) which admits a power series expansion (3.1) Ramanujan's theorem refers then to the formal identity [29] ∞ and implies that the Mellin transform of F (x) is given by The function λ(s), extended over the full complex s-plane, can thus be simply obtained from the discrete n-functional dependence of the λ(−n) coefficients of the Taylor expansion of F (x) by the formal replacement n → −s. The proof of this beautiful theorem was provided by Hardy [30] and it is based on Cauchy's residue theorem as well as on the Mellin-Barnes representation. The basic assumption in Hardy's proof is a growth restriction on |λ(s)| which assures that the series λ(0) − λ(−1)x + λ(−2)x 2 − λ(−3)x 3 + · · · has some radius of convergence. In our case F (x) will be the hadronic photon self-energy function Π(Q 2 ), with x ≡ Q 2 t0 , and Hardy's growth restriction is equivalent to the one required to write a dispersion relation for Π(Q 2 ).
At small Q 2 values, the hadronic photon self-energy function Π(Q 2 ) in QCD has indeed a power series expansion:  A simple example of this procedure was discussed in ref. [2] in the case of vacuum polarization in QED at the one loop level where, in that case, the corresponding Mellin transform is exactly reproduced from its knowledge at just three s values: e.g. s = 1, 0, and −1.

III.2 Marichev's Class of Mellin Transforms.
The class in question is the one defined by standard products of gamma functions of the type with constants C, a i , b k , c j and d l and where the Mellin variable s only appears with a ± coefficient. The interesting thing about this class of functions is that all the Generalized Hypergeometric Functions have Mellin transforms of this type [31]. As a result, many functions have a representation in terms of Mellin-Barnes integrals involving linear combinations of standard products of the Marichev type in Eq. (3.6). 2 In our case, the monotonicity property in Eq. (2.7) of the QCD Mellin transform implies precise restrictions on the subclass of Marichev-like functions that one must consider when trying to implement successive approximations. In that respect we have been particularly helped by some relatively recent mathematical literature [33,34,35]. The authors of these references have studied the general conditions for the convergence of a very general class of Mellin-Barnes integrals, which include those of the Marichev class, and their results can be summarized as follows.
Consider the rather general type of Mellin-Barnes integral . (3.7) In our case this will apply to the Mellin-Barnes integral in Eq. (2.10) where Quite generally, the authors of refs. [33,34] have studied the properties of the mapping which integrals like those in Eq. (3.7) establish between the Mellin s-plane and the z-plane. This is illustrated in Fig. (2) where the crosses denote the positions of the poles in the integrand of Eq. (3.7): in blue the poles at the left of the fundamental strip (represented by the green strip in the figure) and in red at the r.h.s. of the fundamental strip. In the z-plane we show the disc |z| ≤ R in blue, with R the radius of convergence, and the cut starting at Re(z) ≥ R 3 . The converse mapping theorem of ref. [28] relates in a precise way the singularities in the complex s-plane of the integrand in Eq. (3.7) to the asymptotic expansions of I(z) for z large (the red mapping in Fig. (2)) and for z small (the blue mapping in Fig. (2)). Following refs. [33,34,35] we are instructed to consider the two quantities: Then, the region where the integral I(z) converges is | arg z| < π 2 α (see e.g. [33]), and there are three cases to be considered [34,35]: • If ∆ > 0, closing the integration contour to the left leads to a series representation of the integral I(z) which converges for any value of z, but closing the contour to the right gives a divergent asymptotic expansion.
• If ∆ < 0, closing the contour to the right leads to a series representation of I(z) which converges for any value of z, but closing the contour to the left gives a divergent asymptotic expansion.
• If ∆ = 0, closing the contour to the left and to the right gives two convergent series, the first series obtained by closing to the left converges within a disk |z| < R whereas the other one converges outside this disk. Moreover, if α > 0, the two series are the analytic continuation of each other.
These three cases are illustrated in Fig. (3).
We are now in the position of fixing the class of successive Mellin approximants M N (s) that we should use to ensure that they converge in the same way as the full QCD Mellin transform M(s) does.

Figure 3:
Behaviour of the series expansions of I(z) depending on the sign of ∆ for |z| < R (the blue region) and |z| > R. The label div. denotes the regions where the asymptotic expansion is divergent or does not exist. The cut is represented by the green zigzag line.
Associated to each M N (s) approximant there will be a corresponding Π N (Q 2 ) approximant to Π(Q 2 ) (via Eq. (2.10)) and, therefore, a corresponding a HVP µ (N ) approximant to a HVP µ (via Eq. (2.17)). The input will be that we know the values of the first few moments including their errors and their correlation matrix, either from a LQCD determination or from a dedicated experiment. Given this input, we shall then restrict the successive Marichev-like Mellin approximants in Eq. (3.6) to those satisfying the following criteria: In practice, due to the fact that the sequence of poles from Γ(a i − s) is at s = a i + n and the one from Γ(c j + s) at s = −c j − n with n ∈ N implies the restrictions: Re a i ≥ 1 and Re c j ≥ 0 . 3. We also want the corresponding Π N (Q 2 )-function (see Eq. (3.23) below) to the Mellin approximant M N (s) to converge for z ≡ Q 2 t0 both for |z| < 1 and |z| > 1 which, according to the convergence conditions discussed above, requires that (3.14) 4. Finally, we want the two series generated by the Π N (Q 2 ) approximant for |z| < 1 and |z| > 1 to be the analytic continuation of each other which implies This, combined with Eq. (3.14), implies l max < 1 and hence the absence of Γ(d l + s) factors in the denominator of Eq. (3.6).
From the above considerations we conclude that, in the case of HVP in QCD, the only Mellin approximants of the Marichev class that one must consider are those restricted to the subclass: Furthermore, the monotonicity property of the QCD Mellin transform requires that (see e.g. ref. [36]) which implies the asymptotic behaviour each term must satisfy the restrictions in Eqs. (3.17) and (3.18) with real constants C N1 , C N2 , · · · such that Besides the matching to the input moments in Eq. (3.11), all the MBa that we shall use will be constrained to satisfy the leading pQCD short-distance behaviour 4 The Generalized Hypergeometric Function [37] is defined, for |z| < 1, by the series where in the second line we use the Pochhammer symbol This series has P numerator parameters, Q denominator parameters and one variable z. Any of these parameters are real or complex, but the b parameters must not be negative integers. The case where P = 2 and Q = 1 corresponds to the so called Gauss Hypergeometric Function. The sum of this type of series, when it exists, defines a Generalized Hypergeometric Function (GH-Function).
The reason why we are interested in GH-Functions is that, inserting the general expression in Eq.
which is given by the series in Eq. (3.25) for | Q 2 t0 | < 1, with its analytic continuation defined by the underlying Mellin-Barnes integral, Eq. (3.23) in this case. The corresponding Adler function is also a GH-Function: The reason why we are interested in Meijer's G-Functions is that the inverse Mellin transform of M N (s) corresponding to Eq. (2.6), i.e. the Mellin Barnes integrals , (3.31) and have the property that For the class of Marichev-like M N (s) functions in Eq. (3.16) this results in a set of equivalent spectral functions: These successive equivalent spectral functions, alike the physical spectral function, are only defined for t ≥ t 0 but they are not expected to reproduce, locally, the detailed physical shape unless the level of approximation reaches the exact solution (as it is the case in the QED example at the one loop level discussed in ref. [2]). However, when inserted in a dispersion relation integral, they reproduce the predicted smooth behaviour of the successive self-energy functions Π N (Q 2 ) and Adler A N (Q 2 ) functions. It is in this sense that we call them equivalent.
The explicit form of these general expressions for the first N = 1 and N = 2 cases are as follows: • N=1 This corresponds to the case where we only know the first moment M(0). Then to ensure the pQCD pole behaviour at s = 1. The only free parameter b 1 is then fixed by the matching condition M 1 (0) = M(0) and one finds In this simple case the equivalent spectral function is (3.38) • N=2 This corresponds to the case where we know the first two moments M(0) and M(−1). Then and the parameters a 2 and b 2 fixed by the two matching conditions the corresponding Adler function is and the equivalent N = 2 spectral function is 6 : We next propose to show the application of the Mellin-Barnes approximants discussed above to a non trivial example.
We wish to test the techniques developed in the previous section with a more complicated example than the lowest order QED vacuum polarization discussed in ref. [2]. We suggest to examine the case of the QED vacuum polarization at two loops. The proper fourth order QED spectral function was first calculated by Källen and Sabry in 1955 [39] and later on in ref. [40]. It is given by the following expression: With m the lepton mass in the QED VP-loop and  The asymptotic behaviours of this spectral function are Notice that the behaviour at threshold t ∼ 4m 2 is rather different to the one at the one loop level [2] and the shape of the spectral function, which is shown in Fig. (4), is also very different. and shares with QCD the property of being a monotonously increasing function from s = −∞ to s < 1.
The real part of the fourth order vacuum polarization in QED is also known analytically [39]. It is a rather complicated expression and, therefore, it is a good test to see how well it is approximated by the successive GH-Functions in Eq. (3.28). The shape of the Π QED 4th (Q 2 ) function in the Euclidean is shown in Fig. (6).
We shall discuss this 4th order QED example in a way as close as possible to the QCD case which we shall later be confronted with. Therefore, the input will be the successive values of the moments of the spectral function, i.e. of the derivatives of Π QED 4th (Q 2 ) at Q 2 = 0. The first few Mellin moments  We can now proceed to the construction of a successive set of MBa's to M QED 4th (s) of the type shown in Eq. (3.16) and to the evaluation of the corresponding GH-function approximation to Π QED 4th (Q 2 ) of the type shown in Eq. (3.28). At each approximation step we shall then evaluate the corresponding contribution to the anomalous magnetic moment of a fermion of mass m induced by the 4th order vacuum polarization generated by the same fermion (see the corresponding Feynman diagrams in Fig. (7)), and compare it with the exact result which is known analytically [41]: The result in Eq. (4.7) is a rather complicated expression involving higher transcendental numbers with important numerical cancellations among the different terms and, therefore, it should provide a good test. We want to investigate how well we reproduce this exact result using the Mellin-Barnes integral representation in Eq. (2.17) which, when adapted to this case, reads as follows: which must be singular at s = 1. This fixes the a parameter to a = 1 and the overall normalization to This corresponds to the case where we know the slope and curvature of Π QED 4th (Q 2 ) at Q 2 = 0, i.e. M QED 4th (0) and M QED 4th (−1). This information is similar to that already available from LQCD 7 . We shall therefore discuss it in detail.
The Mellin approximant in this case has two parameters a and b: 13) and the leading short-distance constraint fixes the overall normalization to with the parameters a and b fixed by the two matching equations: 15) or equivalently .
gives the result  At this stage it is also interesting to compare the exact Mellin transform shown in Fig (5) with the one corresponding to the N = 2 approximation. This is shown in Fig. (9) where the blue dotted curve is the N = 2 approximation. The agreement of the two curves down to s −3 is quite remarkable. In order to see the difference between these two curves we show in Fig. (10) the plot of their ratio. The M 2 (s)/M(s) ratio turns out to be greater than one everywhere, except in the interval −1 ≤ s ≤ 0. This is why the N = 2 result approaches the exact value of the anomaly from above. The quality of the interpolation between s = 0 and s = −1 provided by the N = 2 approximation is shown at the right in Fig. (10). Notice the scale in the figure, e.g. the value at the minimum of the ratio shown in this figure is 0.9937 compared to one.

Figure 11:
The red curve is the exact 4th order QED VP-function. The dotted blue curve is the N = 2 approximant. Both curves are shown in α π 2 units.
According to Eq. (3.28), the N = 2 GH-function approximant to Π QED 4th (Q 2 ) is given by the expression (z ≡ Q 2 4m 2 ): where 3 F 2 1 1 a 2 b − is the GH-Function defined by the series: and a and b have the values given in Eq. (2.17). Figure (11) shows how well the MBa Π QED (N =2) (Q 2 ) (blue curve) does when compared to the exact function (red curve). From this comparison, one can qualitatively understand why the N = 2 approximation already reproduces the exact value of a VP in Eq. (4.7) at the 0.5% level.
The equivalent spectral function corresponding to the N = 2 approximation is given by the Meijer's G-Function: and its shape, compared to the exact spectral function, is shown in Fig. (12). Notice that the equivalent spectral function corresponding to the unique Padé approximant constructed with M QED 4th (0) and M QED 4th (−1) would be just a delta function.

IV.1.3 The N = 3 MBa.
This corresponds to the Mellin approximant with , (4.27)  Fig. (14). Notice the scale in the left plot of Fig. (14) as compared to the one in Fig. (10) and the improvement in the figure at the right which is plotted at the same scale as Fig. (10). An accuracy of 0.004% is already much beyond what is required of the HVP contribution to the muon anomaly in QCD, but for the sake of testing the approximation procedure that we are advocating, let us try further possible improvements.
, (4.32) and the four parameters a 1 , a 2 , b 1 and b 2 solutions of the matching equations: ,  .16) is an excellent method to approach, rather quickly in this case, the exact result with an excellent accuracy. The question which, however, arises is: how far can one go?. The exact Mellin transform of the QED fourth order spectral function, contrary to the second order one discussed in ref. [2], is expected to be a much more complicated expression than just a simple standard product of the Marichev class in Eq. (3.16). Therefore, a priori, one expects these approximations to break at some N -level where no acceptable solutions exist any longer. Let us then proceed to examine what happens when one tries higher N -approximants of a single standard product.
IV.1.5 The N = 5 MBa. 40) and the parameters a 1 , a 2 , b 1 , b 2 , b 3 solutions of the matching equations: . which reproduces the exact value at the level of 0.00018%, still an improvement with respect to the N = 4 Approximation! This is, however, the best one can do in the two loop QED case with single Mellin approximants of the type shown in Eq. (3.16). Indeed, if one tries to improve with a N = 6 approximant of this type, one finds that all the solutions for the parameters a 1 , a 2 , a 3 , b 1 , b 2 , b 3 from the matching equations bring in complex numbers with real parts which are inside of the fundamental strip, in contradiction with the initial requirements for an acceptable solution that we imposed. This is the signal that, in our example, single Marichev-like approximants break down at a critical N -level where the function Π QED 4th (Q 2 ) cannot be approximated any longer with just one GH-Function. It is possible, however, to extend the class of approximants to superpositions of standard products as indicated in Eq. (3.16) and in fact this is what we shall do in the case of QCD.
From the previous analysis we conclude that, in the case of the QED fourth order vacuum polarization, the best prediction we can make with single Marichev-like MBa's is an average of the N = 4 and N = 5 approximants with an error estimated from the deviation of this average to the N = 4 and N = 5 results i.e., V Test of MBa with experimental HVP Moments.
The KNT collaboration [9] has kindly provided us with the values of the first few moments of the hadronic spectral function with their errors, as well as their covariance matrix. These moments were obtained using the same hadronic spectral function which results in the second number quoted in Eq. (3.16). It provides us with a good test of how well the approximants that we propose work when applied to a set of hadronic moments with realistic errors. The first five moments with their errors are given in Table (2) and their correlation matrix is given in Table (3) in the next section. We observe that the relative errors of the first two moments M(0) and M(−1) in Table (2) are smaller than the relative error in the determination of the lowest order HVP contribution to a HVP µ in Eq. (1.3) [9]. The higher moments M(−n) for n = 2, 3, ... have higher relative errors but they of course contribute less and less to the total a HVP µ determination. We shall next proceed, like in the previous section, to the construction of successive MBa's of the type shown in Eq. (3.16) and to the evaluation of the corresponding GH-Functions Π QCD N (Q 2 ) and 1 π ImΠ N (t). At each approximation we shall then evaluate the corresponding a HVP µ (N ) contribution to the muon anomlay. In the next subsection we shall only consider as input the center values of the moments in Table (2) and postpone the error analysis for later discussion in the next subsection. This corresponds to the MBa which one can construct when only the first moment M(0) is known. In this case where the singularity at s = 1 is the one associated to the asymptotic leading behaviour of the QCD spectral function with u, d, s, c, b and t quarks in Eq. (2.1). Matching the value of M 1 (s) at s = 0 with the one from the experimental determination in Table ( Table (2) with their errors, which are too small to be seen at the scale in the figure. The agreement, at the precision of the scale of the figure, is excellent.
Inserting the expression of the first Mellin approximant M 1 (s) in the integrand at the r.h.s. of Eq. (2.14) gives the result of the first MBa to the muon anomaly: which reproduces the central value result in Eq. (1.3) [9] surprisingly well: to 0.8%. In order to understand why the N = 1 MBa is already so good, let us explore more in detail the plot of M 1 (s) in Fig (15). To better observe the deviations between the experimental moments and the predicted moments we plot in Fig. (16) their ratio as a function of s = −n, n = 0, 1, 2, . . . . The red curve shows the shape of the N = 1 MBa in Eq. (5.1).
The blue circles are the experimental values in Table (2).  Figure 16: Plot of the ratio of the experimental moments in Table (2) with their errors to those predicted by the N = 1 Mellin-Barnes-Approximation.
The deviation of this ratio from one shows the discrepancy. Notice that, here, only the value of the M(0) moment has been used as an input. The predicted values of M(−1), M(−2) and even M(−3) turn out to be rather close to the experimental values, although already the predicted M(−3) and certainly the predicted higher moments are not compatible with the experimental statistical errors. Higher moments, however, contribute less ans less to the total value of the anomaly and this is why a HVP µ (N = 1) turns out to be already such a good approximation. Why does the N = 1 MBa do a better job in the case of QCD than in the two loop QED case we discussed before? The reason for this is that in the QCD case, contrary to the QED case, there are resonances in the low energy region of the spectral function with mass scales which, relative to the muon mass, enhance the contribution of the low moments, in particular M(0). If instead of the muon anomaly we were considering the electron anomaly, the N = 1 MBa would already be giving a result with an accuracy comparable to the full determination.
Although, given the result in Eq. (5.3) and the present accuracy from experiment, there seems to be little room for improvement, let us examine what happens when one tries the N = 2 MBa.
Here the Mellin approximant has the analytic form The shape of the M 2 (s) Mellin transform turns out to be rather similar to the M 1 (s) one in Fig. (16). In order to appreciate the differences between the N = 1 and N = 2 MBa's, we compare in Fig. (17) the ratios of the experimental moments to those of the M 2 (s) prediction (the red dots) and to those of the M 1 (s) prediction (the blue dots). The overall shape of the red dots is clearly better because they are nearer to one.  Figure 17: Plot of the ratio of the experimental moments in Table (2) with their errors to those predicted by the N = 2 MBa in red and the N = 1 MBa in blue.
With the expression of the second Mellin approximant M 2 (s) inserted in the integrand at the r.h.s. of Eq. (2.14) we get as a result of the N = 2 MBa to the muon anomaly: As discussed in the previous section, the MBa technique allows to reconstruct as well Π N (Q 2 ) approximants of the HVP self energy in terms of GH-functions. The corresponding N = 2 approximant is (z = Q 2 t0 ): with a 1 and b 1 given in Eq. (5.7). The shape of the function Π QCD N =2 (Q 2 ) is shown in Fig. (19). Plots of the spectral function associated to the N = 2 MBa are also shown in Figs. (20). Although, asymptotically, the N = 2 MBa spectral function approaches the pQCD value it can only be considered a smooth interpolation of the physical spectral function which, as we know, has a lot of local structure. This interpolation, however, when inserted in the r.h.s. of Eq. (1.4) reproduces the determination of the anomaly using the experimental spectral function at the 0.5% level already mentioned. It is in this sense that it is a good interpolation.
We shall next explore what happens when one tries to improve the N = 2 MBa with higher approximants and further input from the experimental values of higher moments.
The corresponding Mellin approximant which generalizes the one in Eq. (5.1) has the analytic form  (5.13) and the equivalent one with b 1 b 2 . These "solutions", however, are not acceptable because they generate a pole at s = a 1 which is inside of the fundamental strip in contradiction with first principles, as discussed in Section III.3. Nevertheless, the negative numerical values of a 1 and b 2 are in fact rather close to each other. Had they been exactly the same, there would have been a cancellation between Γ(a 1 − s) and Γ(b 2 − s) in Eq. (5.11) indicating that it is not possible to improve beyond N = 2 with a single Marichev-like function. The situation here is rather similar to the one encountered earlier when considering the N = 6 MBa in the QED example.
The fact that in QCD the simple Marichev-like approximants fail to find physical solutions already at the N = 3 level is perhaps not so surprising. One does not expect, beyond a certain level of accuracy, to be able to approximate Π QCD (Q 2 ) at all Q 2 values with just one GH-function. One may, however, ask: is it possible to find generalizations of the simple Marichev-like MBa's which, when using more than the first two moments in Table (2) as an input, provide acceptable solutions to compare with a HVP µ in Eq. (1.3) [9]? As already mentioned at the end of Section IV there is a positive answer to that. It consists in using standard superpositions of Mellin approximants of the type indicated in Eq. (3.16). This, in turn, implies specific superpositions of GH-Functions which approximate the self-energy Π QCD (Q 2 ) in the Euclidean, and hence a HVP µ .
The simplest superposition which gives acceptable solutions to the matching equations, when one knows three moments in the HVP case, consists of the sum of one N = 2 MBa and one N = 1 MBa: Plots of the N = 2 + 1 Adler Function versus z = Q 2 t0 .
The corresponding sum of HG-Functions to the M 2+1 (s) MBa in Eq. (5.15) which results as an approximation to the HVP self-energy is now 18) and the corresponding approximation to the Adler function is The shape of this Adler function is shown in Fig. (21).
With the first four moments of HVP as an input, there is a new superposition of MBa's which gives an acceptable solution to the matching equations. It is the following linear combination of a N = 2 MBa and two N = 1 MBa's: .   Table (2) to those of the N = 2 + 1 + 1 MBa.
Notice the difference of scale in the vertical axis, as compared to the one in Fig. 17.
The Adler function associated to M 2+1+1 (s) in Eq. (5.20) is the sum of three GH-Functions: and its shape is shown in Fig. (24). Plots of the spectral function corresponding to the N = 2 + 1 + 1 MBa are also shown in Fig. (25). The plots already exhibit underlying features of the hadronic structure.  Table (2) and their correlation matrix is given in Table (3). One can see that the values of these moments are highly correlated, reflecting the fact that they all have been extracted from different integrals of the same input data on the spectral function.
The statistical part of the analysis is standard. We first construct the covariance matrix C ij of the first N moments obtained from experiment M(1 − i) , i = 1, . . . , N :  where ρ ij is the correlation coefficient between the moment #i and the moment #j, each with Gaussian uncertainty σ i and σ j . Then we define a χ 2 function associated to a given Mellin-Barnes approximant M N (s), which depends on a set of parameters (a k , b k ): 27) and minimize this χ 2 with respect to the set of parameters (a k , b k ). The errors are sufficiently small to ensure that a point-like estimate is an excellent approximation, and we obtain the covariance matrix in the (a k , b k ) parameter space from the Hessian matrix of the χ 2 function computed at its minimum.
Using linear error propagation we can then calculate the statistical uncertainty on a HVP µ , as reported in the third column of Table (4). The fact that all the approximants have a similar uncertainty that coincides with the one of the complete evaluation of a HVP µ [9] is a sign that the statistical information is saturated by all our MBa's.
Our results would not be complete without a study of the systematic shift associated to the   The pink band is the full experimental result of ref. [9].
successive MBa's which interpolate the values of the experimental moments and reconstruct the full Mellin functions. With this aim, in addition to the MBa's discussed in detail in the previous section, we have also tested alternative parameterizations for N = 2, 3, 4 which are obtained by changing the location of the poles in the superposition terms ( e.g. Γ(2 − s) instead of Γ(1 − s) in Eq. (5.14)). These alternative MBa's have also valid solutions for the corresponding (a k , b k ) parameters and, therefore, can also be considered as good alternative choices. The results of all the evaluations of a HVP µ which we have made are plotted in Fig. (26), as a function of the number of input moments N . We observe that the successive results converge towards the experimental value in Eq. (1.3).

VI Conclusions and Outlook
Equation (2.4) shows that moments of the hadronic spectral function are equivalent to derivatives of the hadronic self-energy function Π(Q 2 ) at Q 2 = 0. The latter are accessible to LQCD simulations as well as to eventual dedicated experiments. We have shown how, from an accurate determination of the first few moments, one could reach an evaluation of the HVP contribution to the muon anomaly with a competitive precision, or even higher, than the present experimental determinations.
The method that we propose uses a new technique of Mellin-Barnes approximants which has been explained and justified in detail in the text. Essentially it is based on generic QCD properties which fix the class of Mellin transforms M(s) of the spectral function that one can use as successive approximants. The muon anomaly a HVP µ , in terms of these M(s)-functions, is given by the Fourier transform in Eq. (2.17). The corresponding approximations to the hadronic self-energy function Π(Q 2 ) are well defined Generalized Hypergeometric Functions which we have given explicitly and the approximations to the spectral function are also given in terms of Meijer's G-Functions. This offers the possibility of applying the same techniques developped here to the case where the information from LQCD, or from experiment, is given in terms of determinations of the self-energy function Π(Q 2 ) at fixed Euclidean Q 2 -values, as e.g. in ref. [15]. We plan to discuss this in the near future.
We have illustrated the practical application of the method with the example of the QED contribution to the muon anomaly from the vacuum polarization Feynman diagrams in Fig. (7). We have also discussed the case where one uses as an input the experimental values of the first moments provided to us by the collaboration of ref. [9]. We find that, in this case, our approach reproduces very well their complete phenomenological analysis.
We shall first show how performing the integral in the r.h.s. for N = 1 reproduces the l.h.s. For that we make a choice of s with Re(s) ∈]0, 1[, e.g. s = 1 2 + iτ . Then LQCD determinations of Mellin Moments should be consistent with these constraints.