A Hybrid Strategy for the Lattice Evaluation of the Leading Order Hadronic Contribution to $(g-2)_\mu$

The leading-order hadronic contribution to the muon anomalous magentic moment, $a_\mu^{\rm LO,HVP}$, can be expressed as an integral over Euclidean $Q^2$ of the vacuum polarization function. We point out that a simple trapezoid-rule numerical integration of the current lattice data is good enough to produce a result with a less-than-$1\%$ error for the contribution from the interval above $Q^2\gtrsim 0.1-0.2\ \mathrm{GeV}^2$. This leaves the interval below this value of $Q^2$ as the one to focus on in the future. In order to achieve an accurate result also in this lower window $Q^2\lesssim 0.1-0.2\ \mathrm{GeV}^2$, we indicate the usefulness of three possible tools. These are: Pad\'{e} Approximants, polynomials in a conformal variable and a NNLO Chiral Perturbation Theory representation supplemented by a $Q^4$ term. The combination of the numerical integration in the upper $Q^2$ interval together with the use of these tools in the lower $Q^2$ interval provides a hybrid strategy which looks promising as a means of reaching the desired goal on the lattice of a sub-percent precision in the hadronic vacuum polarization contribution to the muon anomalous magnetic moment.


Introduction
Current determinations of a µ = (g−2) µ /2 in the Standard Model (SM) show a discrepancy of about 3σ with respect to experiment [1,2]. Moreover, a new Fermilab experiment expects to reduce the error by a factor of 4 in the near future which, should the central values stay the same, would mean a deviation from the SM value of about 5σ! Clearly this problem requires attention. 1

Speaker
At present, the largest component of the error on the SM prediction is that on the leading order hadronic vacuum polarization contribution, a LO,HV P µ . Since the relevant scale for this contribution is m 2 µ , this requires a difficult nonperturbative calculation in QCD. Fortunately, dispersion relations allow a determination of this contribution by relating the vacuum polarization diagram to the e + e − hadroproduction cross-section. However, the discrepancies between different experiments in the most arXiv:1410.8405v1 [hep-lat] 30 Oct 2014 relevant channel, e + e − → π − π + , [3][4][5][6] 2 , together with the need to maximally reduce the error in view of the coming Fermilab experiment has spurred the community to provide an independent determination of a LO,HVP µ from first principles on the lattice [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Such a calculation will have the extra benefit of being a very good testing ground for more difficult problems such as the light-by-light contribution to a µ , for which, at present, there are only model estimates [2].
A convenient representation for the calculation of a LO,HVP µ on the lattice is given by [8,24] a LO,HVP where m µ is the muon mass and is the subtracted polarization, defined from the hadronic electromagnetic current twopoint function, Π µν (Q), via The value of a LO,HVP The two-point function Π µν (Q) can in principle be computed on the lattice for non-zero Q whence Π(Q 2 ) can be extracted. However, this extraction is complicated in practice. A typical situation is depicted in Fig.  1. In this figure, the blue dashed curve shows the result of a physically motivated model for the vacuum polarization in the I = 1 channel based on the experimental data obtained in non-strange hadronic τ decays, supplemented by a successful model of duality violations at energies above the τ mass [25][26][27][28][29][30] (see Ref. [31] for more details). Also plotted is fake lattice data at a set of Q 2 i values corresponding to a recent lattice simulation on a 64 3 × 144 lattice, with a = 0.06 f m, m π = 220 MeV and periodic boundary conditions [32]. The fake data is obtained by letting the model Π(Q 2 i ) fluctuate according to the covariance matrix of this simulation . So, Fig. 1 is a realistic representative picture of the situation on the lattice. 2 A useful overview of the experimental situation is given in Figs. 48 and 50 of Ref. [5]. (1) as a function of Q 2 in a τ-data-based model. The peak is at Q 2 ∼ m 2 µ /4 0.003 GeV 2 . Red data points: a set of typical lattice data.
Obviously a direct numerical integration of the data shown in Fig. 1 is not an option for an accurate evaluation of the area under the curve needed in Eq. (1). It is necessary to use some functional form in a fit and make sure that this functional form will faithfully reproduce the curve in Fig. 1, and this is why having a model becomes very useful: by knowing the answer (in the model) we can assess the size of the systematic error made. The figure of merit to keep in mind is that we need a determination of a LO,HVP µ with better than 1% precision.
One technical point: experimental spectral data can only determine the vacuum polarization through a subtracted dispersion relation. So, we can think of our model as one in which Π I=1 (0) = 0. However, to really mimic the situation on the lattice, where the value of Π I=1 (0) is unknown, we will always consider Π I=1 (0) as a free parameter to be determined by fits to the data. The extent to which this value deviates from 0 then quantifies the systematic uncertainty in the determination of Π I=1 (0).

Hybrid strategy: "Divide and Conquer"
Due to the prominent peak at low Q 2 seen in Fig. 1, the integral in Eq. (1) is largely dominated by the contribution at low energies. This can be seen in Fig. 2, where more than 90% of a LO,HVP µ is accumulated below Q 2 ∼ 0.2 GeV 2 . In fact, it is only in this lowenergy region that lattice data is so scarce that fitting is required. For higher energies a naive numerical integration of the data using the trapezoid rule is sufficient [19]. This can be seen in Fig. 3 where, for each value of Q 2 min , we show the systematic uncertainty in a LO,HVP for Q 2 min 0.1 − 0.2 GeV 2 , both uncertainties are well below 1%. Of course, we only know the systematic uncertainty because we know the exact result in the model. This assumes that the subtraction constant Π(0) is known with sufficient precision. In a previous work [19], we showed that a fit of a simple [1, 1] Padé 3 to the fake data in the interval between 0 and 1 GeV 2 is able to determine Π I=1 (0) with an uncertainty, δΠ I=1 (0), smaller than 0.001. An uncertainty δΠ I=1 (0) produces a corresponding uncertainty on the contribution to a LO,HVP Fig. 4 shows this uncertainty. We see that the error remains safely below 1% for Q 2 0.1 − 0.2 GeV 2 . This error will have to be carefully monitored in the final analysis, however, due to the rapid increase at low Q 2 seen in Fig. 4.
These observations tell us that a hybrid strategy in which one divides the integration interval into two parts could achieve the desired precision: the first one covering 0 ≤ Q 2 ≤ 0.1 − 0.2 GeV 2 , and the second one covering 0.1 − 0.2 GeV 2 ≤ Q 2 < ∞. It is only in the first part that one should fit to a functional form, while in the second part one may use a trapezoid rule integration of the data. There is no need for the long (and dangerous) extrapolation down to Q 2 ∼ m 2 µ from the region of "good data" at Q 2 ∼ 1 GeV 2 , as has been customarily done until now. 3 See below for further discussion on Padés.

The low-Q 2 region
In this region of Q 2 we would like to propose a strategy based on three independent tools: Padé Approximants, polynomials in a conformal variable and an NNLO Chiral Perturbation Theory (ChPT) representation supplemented by an analytic Q 4 term.
Padé Approximants are ratios of polynomials whose coefficients are matched onto an equal number of derivatives of the original function Π(Q 2 ) at a single point, usually at Q 2 = 0 [33], or at different values of Q 2 [34]. In the first case they are called one-point Padés, while in the second case they are multi-point Padés. Because the vacuum polarization is a so-called Stieltjes function, there are convergence theorems which control the approximation of Padés to the original function everywhere in the Q 2 complex plane, except right on the cut of the vacuum polarization. For example, for Q 2 > 0, which is the region of interest in the integral of Eq. (1), one has for one-point Padés that where [M, N] H represents the ratio of a polynomial of degree M over a polynomial of degree N, matched onto M + N coefficients of the Taylor expansion of Π I=1 (Q 2 ) about Q 2 = 0. The subscript "H" refers to the fact that these Padés actually denote −Q 2 times the Padés conventionally used in mathematics [33,34]. In Ref. [21] it was shown how one can determine the Taylor coefficients of Π I=1 (Q 2 ) at Q 2 = 0 from time moments of the Euclidean two-point function of the vector current. To get the n-th term of the Taylor expansion, one needs an accurate determination of the 2n + 2-nd time moment of the two-point function. In Ref. [21] four Taylor coefficients were determined for s, c quarks and then Padés were constructed to approximate the vacuum polarization function for all Q 2 and, consequently, the integral in Eq. (1). The contribution for u and d quarks has not yet been done, and is expected to be a lot harder.
However, as we have seen, it is only necessary for the Padés to approximate the vacuum polarization function in the low-Q 2 window 0 ≤ Q 2 ≤ 0.1 − 0.2 GeV 2 . This is a big advantage as the more restricted Q 2 range means fewer derivatives of the correlator are required for an accurate Padé representation in this region. Our tau-data-based model now allows us to investigate how many of these derivatives are needed to reach a given accuracy for a LO,HVP , as a function of Q 2 max . This is shown in Fig. 5 As an alternative to the Padés constructed from the Taylor expansion at Q 2 = 0, there are also the multipoint Padés [13], whose coefficients are fixed by the  value of the original function (and, if available, also its derivatives) at a discrete set of values in a Q 2 interval. 4 . This means in practice that one has to make a fit. As an example, we have done the exercise of fitting the τ-based model data at the set of points Q 2 = 0.10, 0.11, · · · , 0.20 GeV 2 using the [2, 1] H Padé form, withΠ I=1 (0) as a free parameter. We emphasize that this is not the fake data mentioned earlier, which was based on a real set of lattice data. With periodic boundary conditions, the current lattice data shows too large errors and a too small number of values of Q 2 for this type of fit to be successful. So, even though the present exercise cannot be considered as fully realistic now, it may become feasible in the future thanks to error-reduction techniques [35,36] and new theoretical ideas [12,[14][15][16][17]37]. Even so, we find that it is necessary to go to the [2, 1] H Padé to get down to the sub-percent level in the systematic error of the integral in Eq. (1) from 0 to 0.1-0.2 GeV 2 . The rule of thumb is that, in order to achieve the same level of accuracy, Padés constructed from fitting in an interval require one order more than those obtained by the Taylor coefficients at the origin. See Ref. [31] for more details.
Another possibility is a conformal expansion of the subtracted vacuum polarization. The Taylor expansion ofΠ I=1 (Q 2 ) converges only for |Q 2 | < 4m 2 π , which is too small to be useful. The convergence properties can be improved by going to the conformal variable and then expanding in w. The new functionΠ I=1 (Q 2 (w)) should converge faster because the whole Q 2 complex plane is mapped onto the unit disc, with the Q 2 cut at the boundary. The region of convergence now includes those w corresponding to the whole positive real Q 2 axis, which is what is needed in Eq. (1). As in the case of the Padé Approximants, the conformal expansion can be constructed from the lowest Taylor coefficients of the vacuum polarization function which, in turn, can be determined by the Euclidean time moments of the two-point correlator. With up to the fourth order Taylor coefficient, one can construct a linear, quadratic, cubic and quartic polynomial approximations toΠ I=1 (Q 2 ) and compare with the exact result. The result of this comparison is shown in Fig. 6.
This figure shows that the linear polynomial is clearly insufficient. The quadratic version yields much better estimates forâ LO,HVP µ [0, Q 2 max ], i.e. 0.6% and 1% below the exact model values for Q 2 max = 0.1 and 0.2 GeV 2 , respectively. In the case of the cubic representation, the corresponding errors are 0.02% and 0.04%. These numbers are to be compared to 0.3% and 0.5% for the [1,1] H Padé (which has the same number of parameters as the quadratic polynomial), and 0.06% and 0.2% for the [2,1] H Padé (which has same number of parameters as the cubic polynomial).
Also, as in the case of the Padés, one may consider fitting these conformal polynomials in a interval of Q 2 values to avoid the difficulties of the Euclidean time moment calculation. We have repeated the same exercise we carried out for the Padés and the τ-data-based model above in the interval Q 2 = 0.10, 0.11, · · · , 0.20 GeV 2 and find that, again, polynomials in the conformal variable obtained from fitting require one order more than those obtained from the Taylor coefficients at Q 2 = 0 in order to reach the same accuracy forâ LO,HVP Finally, we would like to comment on a strategy based on ChPT [38,39]. Since the region of interest is given by low values of Q 2 , 0 ≤ Q 2 ≤ 0.1 − 0.2 GeV 2 , in principle ChPT should yield a good representation. However, the highest order of ChPT available is NNLO [40,41] and, unfortunately, this turns out to be insufficient. This is not very surprising: the old success of phenomenological descriptions like Vector Meson Dominance suggest that vector resonances have to play a very important role. However, NNLO is just the first order at which vector mesons appear in the subtracted vacuum polarization function through the associated low-energy constants (LECs). Estimating the size of ρ-induced corrections beyond NNLO, one finds NNNLO corrections, for which no complete calculation exists, are likely to become important already for Q 2 ∼ 0.1 GeV 2 .
Given this state of affairs, as an exploratory exercise, we have supplemented the NNLO result with an extra Q 4 term to construct an incomplete NNNLO form which we call NN LO. The idea is that we want to know whether, with sufficiently good low-Q 2 data, at least this form is capable of representingΠ I=1 (Q 2 ) accurately enough to produce a sub-percent result for a LO,HVP µ [0, 0.1 − 0.2 GeV 2 ]. Using again our τ-data based model we can determine the unknown LECs needed to construct our NN LO representation from the first and second derivatives of the vacuum polarization function in this model, and then compare to the exact result forΠ I=1 (Q 2 ) (see Ref. [31] for more details). The result is shown in Fig.  7. As one can see, the NNLO dotted curve fails to give an accurate representation of the red solid curve, representing the exact result, in the region of interest.
For Q 2 max = 0.1 GeV 2 , we find that our NN LO ChPT value forâ LO,HVP µ [0, Q 2 max ] is 0.6% below the exact value, while for Q 2 max = 0.2 GeV 2 , it is 1.4% below. Although the value at Q 2 max = 0.1 GeV 2 is acceptable, this is clearly worse than the approximation obtained using a [1,1] H Padé, which was also determined from the first and second derivatives at Q 2 = 0. For comparison, NNLO ChPT, which corresponds to removing the extra Q 4 term, yields values 4% and 18% above the exact value, at Q 2 max = 0.1 and 0.2 GeV 2 , respectively. This is clearly insufficiently accurate.
As before, we can also attempt to construct our NN LO function from a fit to a set of Q 2 values in an interval, instead of from derivatives at Q 2 = 0. However, Fig. 7 shows that, unlike the case of Padés, for which a fit in the interval [0.1, 0.2] GeV 2 was in principle possible, in the case of the NN LO representation the fit window would have to be below Q 2 ∼ 0.1 GeV 2 . Accurate lattice data for a dense set of such low values of Q 2 may be harder to get. We conclude that a ChPTbased approach may be useful for some sort of consistency check but it is unlikely to be as useful as Padés or the polynomials in the conformal variable we have described.

Errors for the hybrid strategy and conclusions
The contribution from the vacuum polarization function to the muon a LO,HVP µ , Eq. (1), requires knowledge of this function for all values of Q 2 . Even though for Q 2 large enough, say Q 2 2 GeV 2 , one may apply perturbation theory, the largest contribution comes from values of Q 2 ∼ m 2 µ /4 which are much smaller. Lattice data provides a discrete set of points in this region of small Q 2 but, currently, neither the accuracy nor the density of these points is sufficient to allow a numerical integration to cover the region of the integral going from Q 2 ∼ 2 GeV 2 down to Q 2 = 0. The standard method of calculation until now has been to use an extrapolation of the data from Q 2 ∼ 1 − 2 GeV 2 , where the data is quite accurate, all the way down to Q 2 = 0. This results in an unwanted systematic error which makes a reliable sub-1% precision in the total contribution an impossible goal to reach.
In Ref. [31] we have pointed out that it is advantageous to divide the integration region into two parts, one covering 0.1 − 0.2 GeV 2 Q 2 2 GeV 2 and another one covering 0 Q 2 0.1 − 0.2 GeV 2 . The contributions from these two regions can then be evaluated using a hybrid strategy. First, for the upper part, one can apply a simple trapezoid-rule numerical integration. Existing lattice data is already good enough to produce a value for this part of the integral which is sufficiently accurate. The problematic region is the lower part, where current lattice data shows large errors and the Q 2 coverage is too sparse to allow a numerical integration.
We have pointed out that there are three methods likely to have some utility in dealing with the region 0 Q 2 0.1 − 0.2 GeV 2 . These are Padés, polynomials in the conformal variable and ChPT. Of these, we have seen that Padés and the polynomials in the conformal variable will probably be the most efficient, while ChPT, because it is only fully known up to NNLO, will not be so optimal.
The Pades and polynomials in the conformal variable can both be constructed from the values of the derivatives of the subtracted polarization with respect to Q 2 at Q 2 = 0, values which can, in principle, be determined from the Euclidean time moments of the vector current correlator, as Ref. [21] has shown for s and c quarks. A question that still remains is whether such a determination will be precise enough for the case of u and d quarks. Alternatively, one could also construct these Padés and polynomials in the conformal variable by fitting lattice data in a subset of points in the region 0 Q 2 0.1 − 0.2 GeV 2 , once the data in this region becomes better in the future.
In order to understand how precise the values of Π I=1 (Q 2 ) and its derivatives at Q 2 = 0 have to be to reach the desired sub-percent total error, we can construct the Padé [1, 1] H (which was found to be sufficient to reach this precision in the low-Q 2 region) aŝ Π(Q 2 ) = Π(Q 2 ) − Π(0) = a 1 Q 2 1 + b 1 Q 2 .
Errors δa 1 and δb 1 on the parameters a 1 and b 1 produce associated errors on a LO,HVP µ [0, Q 2 min ]. Using our τ-data-based model we have found that a sub-percent error on a 1 , together with at most a few percent error on b 1 , will be enough to obtain a sub-percent error on a LO,HVP µ [0, 0.1−0.2 GeV 2 ]. Further quantitative studies using our τ-based model will become possible once improved lattice data becomes available. This will allow us to construct fake data sets with realistic errors and correlations from the point of view of the lattice, which can then be used to assess the systematic error associated with the use of fit forms in the low Q 2 region of any lattice evaluation of the contribution from the hadronic vacuum polarization to the muon anomalous magnetic moment.
MG is supported in part by the US Department of Energy, KM is supported by the Natural Sciences and Engineering Research Council of Canada, and SP is supported by CICYTFEDER-FPA2011-25948, 2014 SGR 1450, and the Spanish Consolider-Ingenio 2010 Program CPAN (CSD2007-00042).