Towards a precise lattice determination of the leading hadronic contribution to (g-2)_mu

We report on our computation of the leading hadronic contribution to the anomalous magnetic moment of the muon using two dynamical flavours of non-perturbatively O(a) improved Wilson fermions. The strange quark is introduced in the quenched approximation. Partially twisted boundary conditions are applied to improve the momentum resolution in the relevant integral. Our results, obtained at three different values of the lattice spacing, allow for a preliminary study of discretization effects. We explore a wide range of lattice volumes, namely 2 fm<L<3 fm, with pion masses from 600 to 280 MeV and discuss different chiral extrapolations to the physical point. We observe a non-trivial dependence of a_mu(HLO) on m_pi especially for small pion masses. The final result, a_mu(HLO)=618(64)*10^(-10), is obtained by considering only the quark connected contribution to the vacuum polarization. We present a detailed analysis of systematic errors and discuss how they can be reduced in future simulations.


INTRODUCTION
The magnetic moment of a charged lepton is extracted from the vertex function describing the interaction between the lepton and a photon in the limit of vanishing photon momentum. The corresponding anomalous magnetic moment a l is then defined as half the difference between the gyromagnetic factor g and its classical value of 2, i. e. a l = (g l − 2)/2. In the case of the electron the quantity is dominated by QED contributions. The one-loop result was obtained by Schwinger more than sixty years ago [1], and since then a e has reached an accuracy better than one part per billion on both the theoretical and the experimental sides, which yield results in beautiful agreement (see section 3 of [2], and references therein). The anomalous magnetic moment mediates helicity flip transitions [3], which implies that quantum corrections due to heavier particles of mass M , in the Standard Model or beyond, are proportional to m 2 l /M 2 . For this reason the muon anomalous magnetic moment a µ is regarded as a sensitive probe for effects of nearby New Physics. However, by the same argument, given that m µ ≤ m π , the hadronic contributions to a µ are large and notoriously difficult to quantify. While the experimental and theoretical estimates have each reached similar levels of precision of 0.5 ppm, a tension by 2 or 3 standard deviations between theory and experiment persists [4,5]. Before invoking New Physics as the reason for this tension the theoretical result and, in particular, all contributions due to hadronic effects, must be corroborated. The uncertainty is dominated mainly by the hadronic leading order (a HLO µ ) and secondly by the hadronic light-by-light contributions. Currently a HLO µ is estimated via a phenomenological approach based on the evaluation of a dispersion integral. In the low-energy regime the spectral function in the integrand must be determined experimentally, either from the cross section e + e − → hadrons or from the rate of hadronic τ -decays. Both methods suffer from different systematics [2,5] and yield results in slight tension among each other.
None of them however reduces the discrepancy between theory and experiment on a µ below the two standard deviation level. A purely theoretical estimate of a HLO µ from a first-principles approach is clearly desirable, and the present work represents our first step in that direction by using the lattice regularization of QCD.
The hadronic vacuum polarization contribution to a µ has received considerable attention from the lattice community in recent years. Initial studies have been performed in the quenched approximation [6,7] and in the theory with two [8] and three dynamical flavours [9,10]. Compared to the determination of "standard" quantities such as hadron masses, quark masses and decay constants (see [11] for a review), lattice calculations of a HLO µ are extremely challenging. The relevant lattice quantity (the hadronic vacuum polarization discussed in the next section) receives contributions from quark disconnected diagrams, which are intrinsically noisy and difficult to estimate with good statistical accuracy at a reasonable numerical cost. Moreover, the dependence of the hadronic vacuum polarization on the momentum transfer must be accurately traced down to momenta of order m 2 µ and beyond. This value is well below the lowest Fourier momentum (2π/L) 2 which can be reached in current lattice QCD simulations. As a consequence finite size effects may conceivably be large on results obtained within the "standard" approach. In addition, the vacuum polarization receives sizeable contributions from the low-lying vector resonances. Those should be properly accounted for in simulations at sufficiently light quark masses for the lowest vector meson to be a resonance and including the full dynamics of the strange and charm quarks. Finally, as the contribution to a HLO µ from isospin breaking effects may be of the same size of its present uncertainty, those will have to be included in lattice computations (possibly along the lines discussed in [12][13][14]) for them to have a crucial impact on (g − 2) µ phenomenology.
In [15,16] we have shown how (partially) twisted boundary conditions [17][18][19] can be used to improve the momentum resolution in the connected part of the hadronic vacuum polarization, and we have obtained an estimate of the disconnected contribution in Chiral Perturbation Theory, thus addressing the first two systematic effects discussed above. Here we report on numerical results obtained in that setup and show how the fits to the momentum dependence of the vacuum polarization function get substantially stabilized in this way, thereby improving the accuracy of our estimates. Partial results have already appeared in [20,21]. In addition we present a thorough investigation of all sources of systematic errors, including the modelling of the q 2 -dependence of the vacuum polarization, chiral extrapolations, lattice artifacts and finite volume effects.

DEFINITIONS AND LATTICE SETUP
The Euclidean hadronic vacuum polarization (VP) tensor is defined as For space-like momenta, the relation between Π (N f ) µν (q 2 ) and the lowest order hadronic contribution a HLO µ to the anomalous magnetic moment of the muon has been derived in [6,7,22] and reads (suppressing the index N f )

Lattice regularization
We perform our computation on a subset of the gauge configurations generated within the CLS initiative [23] for two flavours of non-perturbatively O(a) improved Wilson fermions [24] and using the standard Wilson plaquette gauge action. The simulation parameters are collected in Table 1, including the values of the twist angle θ, whose rôle will be explained below. Following [7], we have implemented the one-point-split conserved vector current where U µ ∈ SU (3) represents the gauge link in the positive µ direction, and f is again a flavour index. The lattice version of the vacuum polarization tensor then reads By Fourier transforming the expression above into momentum space one gets where q µ = 2πn µ /L andq µ = 2 a sin aqµ 2 with n µ ∈ 0, 1, . . . L/a − 1. We restrict our attention to the case µ = ν to avoid mixings of the composite field V µ (x)V ν (0) with lower dimensional ones in the limit x → 0. It then follows that Π (N f ) µν (q) fulfils the Ward identitieŝ  [25] and have been determined using the definitions and some of the results in [26][27][28]. The values of the lattice scale must be regarded as preliminary. The masses of the lightest pseudoscalar state m π are taken from [29]. The twist angle θ is applied in the spatial x-direction only.
Consequently on the lattice we can relate the scalar vacuum polarization to the tensor one, by mimicking the continuum relation The scalar vacuum polarization extracted in this way approaches its continuum counterpart with a rate proportional to a, where the O(a) effects appear due to off-shell contributions in Π After performing the Wick contractions in Π µν (x) one realizes that connected as well as disconnected quark diagrams contribute. We neglect here the disconnected terms, as done in most of the existing lattice computations in the literature. In [8] such contributions were included for almost half of the ensembles used and found to be negligible within the quoted errors. Following the arguments in [16], the disconnected diagrams are expected to decrease the value ofΠ (N f ) µν (q) by about 10%, which demonstrates the limited accuracy of lattice calculations without their proper inclusion.

Twisting the connected part of the vacuum polarization
By modifying the spatial boundary conditions on the quark fields entering the vector current it is possible to improve the momentum resolution in Π (N f ) (q 2 ) and to access momenta different from the integer multiples of 2π/L. Namely, imposing the condition is equivalent to boosting the momenta in the quark propagator by θ k /L. This technique can be used to modify and refine the lattice dispersion relation of, for example, charged pseudoscalar mesons, or of any momentum dependent quantity, as long as there are no strong final state interactions [18]. This remains true, up to exponentially suppressed finite size effects, for partially twisted boundary conditions [17][18][19], where only (some of) the valence quarks satisfy twisted boundary conditions, while the sea quarks fulfil periodic boundary conditions. Note, however, that the effect of partial twisting obviously vanishes whenever it is applied to flavour-singlet quantities, as the vector current discussed above.
We have shown in [16] that, by introducing a sufficiently large number of valence quarks degenerate with the u, d . . . flavours, the quark connected and disconnected parts of the correlator Π µν (x) can be rewritten as independent correlation functions in Partially Quenched QCD. In such unphysical theories each quark diagram has an unambiguous field-theoretic expression and well-defined continuum and infinite volume limits. The formulation naturally turns the connected contribution into a correlator of flavour-off-diagonal vector currents and thus twisting can be applied to induce arbitrary momentum. This amounts to a simple modification of Eq. 2.7 (restricted to its connected part) in which q k → q k − θ k /L andq changes accordingly.
In practice we have applied twisting in the x-direction only and remained with periodic boundary conditions in the other directions. The twist angles in Table 1 have been chosen in order to achieve a fixed and large density of q 2 values between consecutive Fourier modes and in order to reach a lowest non-zero q 2 around m 2 µ . For the Fourier modes we have considered all the integers values between (0, 0, 0, 0) and (2, 2, 2, 2) in units of 2π/L.

RESULTS
In order to compute a HLO µ we must obtain a continuous description of Π (N f ) (q 2 ) to determine Π (N f ) (0) andΠ (N f ) (q 2 ) entering the integral in Eq. 2.3. At each value of the lattice spacing a we introduce a maximum momentumq 2 max and adopt several different fit ansätze for the loŵ q 2 region belowq 2 max . We additionally impose a matching condition with perturbation theory atq 2 =q 2 max and hence use the perturbative expression to describe the highq 2 region. The results for a HLO µ at different values of the mass of the lowest pseudoscalar state m π are then extrapolated to the physical point using a functional form inspired by Chiral Perturbation Theory.

Fitting procedure
In order to check for the stability of our results against variations in the fitting procedure we use three different ansätze to describe Π (N f ) (q 2 ) in the region 0 ≤q 2 ≤q 2 max : A) a model-independent (2, 3) Padé ansatz defined by a degree 2 over a degree 3 polynomial inq 2 (n coeff = 6), i.e. , C) a vector dominance-model with two poles 1 . This form (n coeff = 5) has been suggested and used in [10] and reads subject to the additional constraint c 2 < e 2 .
We perform correlated fits. For these to be reliable and to produce a meaningful correlated χ 2 large statistics is needed. In [30] a toy model was used to provide the thumb-rule N > ∼ D 2 , where N is the number of measurements and D the number of degrees of freedom in the fit. We therefore randomly select about 25 values ofq 2 and indeed observe stable correlated fits as well as the absence of numerical problems when inverting the average covariance matrix. The errors on the fit parameters are estimated and propagated to a HLO µ by performing the fits and the integral for each bootstrap sample. For the fit ansatz C we compare in Fig. 1 the value of the smaller mass parameter obtained from the fit with the mass measured from the exponential decay of the vector two-point function. We call the latter naive vector mass, as it assumes the lightest vector meson to be stable at all values of m π considered. While this may not be the case for our most chiral points, for the sake of this qualitative comparison such a definition seems nonetheless adequate. A separate study will have to be devoted to the measurements of the vector resonance parameters following [31,32]. The figure shows overall reasonable agreement, with discrepancies which stretch to 15% at most.
As a constraint on the high-q 2 region we always include perturbation theory in the MS scheme at the scale 2 GeV, using the leading-order expression in α s for Π (N f ) (q 2 ) from [33]. Although the NLO in α s is available from the same reference, for our purposes it is enough to include the leading-order, where in addition no internal fermionic loops appear and therefore it is rather simple to apply the formulae to the case of two light dynamical quarks and a quenched strange. In addition, as it will become clear in the following, the perturbative contribution to a HLO µ estimated in this way is extremely small. To evaluate the perturbative formula we use the non-perturbative Λ MS parameter for N f = 2 from [34] and the non-perturbative renormalization factors in [35][36][37] to relate the lattice quark masses to their values in the MS scheme at 2 GeV. As discussed in [7] the function Π(q 2 ) computed on the lattice and in the continuum dimensional regularization can differ by an un-physical and non-universal integration constant. We fix the value of this scheme dependent constant by matching the perturbative expression for Π(q 2 ) to our non-perturbative data atq 2 =q 2 max . Moreover the perturbative prediction diverges as q 2 → 0 and therefore even after this step perturbation theory can be expected to provide a good description of our data only at sufficiently large values ofq 2 . We require the resulting function not only to be smooth but also to have a smooth first derivative on the whole real axis, and in particular atq 2 =q 2 max . Hence, the matching produces one non-trivial relation among the parameters in the fit ansätze above, reducing their number n par from n par = n coeff to n coeff − 1.
The fit ansatz B always produces reduced χ 2 larger than 5 and we therefore discard it. We monitor the stability of the fit results for a HLO µ in the way shown in the left panel of Fig. 2 where, as an example, the values obtained from fits of type A and C to the measurements performed on the F6 ensemble are plotted. A cut χ 2 /dof ≤ 2.5 is applied in producing the figure. The two functions give compatible results forq 2 max > ∼ 1.5 GeV 2 and in the end we opt for the results from ansatz A on all the ensembles, as this yields in most cases the smallest χ 2 values and the smallest errors on a HLO µ . The values ofq 2 max are chosen differently for different ensembles as theq 2 -ranges explored depend on the lattice spacing and on L/a. Specifically the maximum momentum lies in the range 2.4 GeV 2 < ∼q 2 max < ∼ 5.5 GeV 2 . In order to quantify the improvement brought by the use of partially twisted boundary conditions we repeat in Fig. 2, right panel, the comparison between fit ansätze A and C on the F6 ensemble but restricting the data to the Fourier modes only. Twisting clearly helps in stabilizing the fits and the plateau values are reached "earlier" inq 2 max . This increases our confidence in the fitting adopted procedure. Statistical errors are also slightly smaller in the left panel of Fig. 2 with respect to those in the right one.  We summarize our results on a HLO µ for the two-flavour theory and the case with an additional quenched strange quark in Tables 2 and 3  How statistical errors of order 5-10% come about can be easily understood from Fig. 3. In the left panel we show the result forΠ (2) (q 2 ) together with the fit function (ansatz A) and the perturbative curve for the F6 ensemble. The horizontal axis is divided into four regions. Region I, is defined by 0 ≤q 2 < m 2 µ (not visible on the left panel) where theq 2dependence is only constrained by smoothness requirements on the fit function. There is no direct measurement of the vacuum polarization function here, although we could have tuned the θ-angle in order to penetrate this region. However, the error on the corresponding value of the integrand in Eq. 2.3 would have turned out to be larger than 100%. Region II, where m 2 µ ≤q 2 < 0.2 GeV 2 , is accessible only thanks to the use of partially twisted pert (q 2 ) Π relative contributions from the different regions to the integral in Eq. 2.3. Since region I contributes about 25% and is not constrained by any direct measurement of Π (2) (q 2 ), its uncertainty dominates the overall statistical error on a HLO µ , with an ambiguity which, within this approach, inevitably amounts to about 5%. The uncertainty on the contribution from region II is of the same order but slightly smaller, whereas the uncertainties on regions III and IV are sub-dominant.

Chiral extrapolations
We show our results for a HLO µ in the N f = 2 and the N f = 2 plus a quenched strange quark theories as functions of m 2 π in Fig. 4. The blue curves represent our preferred chiral fits and we will discuss in the following how those have been obtained. We start from the observation that in the mass-region explored a curvature is clearly visible in our data. Indeed a linear fit in that region would fail in producing acceptable χ 2 /dof values. We take the observed curvature into account in three different ways.
• We perform a chiral fit inspired by the functional forms derived at NLO in χPT for different well known chiral corrections. These include a chiral logarithm, which should account for the curvature.
• We re-analyze our data following the alternative procedure presented in [8], which is expected to remove the curvature due to the lowest lying vector resonance.
• We impose a cut on the pion mass (< 400 MeV) and linearly extrapolate our data at the four lightest values of m π in the N f = 2 theory. In doing so we negelect cutoff effects and fit the β = 5.2 and β = 5.3 results simultaneously. Such a procedure produces a value for a HLO µ at the physical point which is well consistent with the result from the χPT-inspired fit ansatz in Eq. 3.5 although slightly lower (namely a HLO µ = 508 (62)).
We do not further discuss the third analysis and present in detail the approaches used for the first two fitting procedures.
Since little is known about the functional form describing the dependence of a HLO µ on the light quark masses, we conservatively adopt a form inspired by χPT to extrapolate our data to the physical point. Namely we use the fit function a HLO µ (m π ) = a HLO µ (0) + B (am π ) 2 + C (am π ) 2 ln((am π ) 2 ) , (3.4) with a HLO µ , B and C as free parameters. Such a functional form can be derived for the connected part of a HLO µ from the expressions for the (connected) vacuum-polarization obtained in [16] at NLO in χPT. In addition, since our current set of ensembles covers a wide range of pion masses at only one value of the lattice spacing (a = 0.063 fm), and in order to avoid mixing of cutoff effects with chiral effects, we decide to extrapolate the β = 5.3 data only. The other data sets in Table 1 are used to asses lattice artifacts and finite size effects. Eventually, extrapolations to the continuum limit at fixed pion masses will have to be performed before the chiral extrapolation.
The resulting curves are shown in Fig. 4. We find that the function describes the whole set of data points quite well, even those which are not included in the fit. In the two-flavour case we obtain at the physical point  which is well consistent with the recent N f = 2 result in [8], whereas the inclusion of a quenched strange quark gives where the errors are purely statistical.
In [8] an alternative extrapolation procedure was proposed, with the aim of reducing the m π -dependence. It can be motivated starting from a vector dominance description of where the m π -dependence is explicitly indicated. The quantity g V is related to the vector decay constant f V by g V = f V /m V . The dependence of g V on m π is neglected based on the numerical observation in Fig. 1 of [39]. Eq. 3.7 then leads to .

(3.8)
It is easy to see that, by rescaling the argument of the function f from q 2 to hq 2 , one obtains .

(3.9)
Under these assumptions the choice h = m 2 ρ m 2 V (mπ) in [8] would therefore remove the dependence on m π (up to the mild one in g V ) producing in addition the right physical result since h → 1 in that limit. Such a rescaling has the further advantage of making the dimensionless quantity a HLO µ,h independent from the lattice spacing [8], as opposite to a HLO µ where the lattice spacing value is needed to convert m µ in the weight-function f (see Eq. 2.4) to lattice units.
In order to estimate the uncertainty in our chiral extrapolation we perform the analysis described above, using the naive vector mass to define h. Results for a HLO µ,h are shown in Fig. 5. We expect from the observation in the previous section on the high χ 2 /dof values produced by the fit ansatz B (single vector dominance) that some dependence on m π will remain. However, the rescaling clearly renders such a dependence mild, its effect being particularly strong for heavy pion masses. The curvature visible in Fig. 4 has almost completely disappeared and we therefore linearly extrapolate the results from the three smallest pion masses at β = 5.3 to the physical point. There is no obvious theoretical reason why a non-linear dependence could not survive, however, we do not have a large sensitivity to those terms in the m π -range explored here. We obtain in this way the following results at the physical point which are consistent with those in Eqs. 3.5 and 3.6. Notice that, as a consequence of the rather moderate chiral extrapolation performed, they are also compatible with the values of a HLO µ,h directly estimated at our most chiral point (F7 ensemble).

Residual cutoff and finite size effects
The set of our simulations covers a rather wide range of lattice volumes and lattice spacings.
We are therefore able to address the issue of finite volume effects and cutoff effects. This study is not yet complete and should be extended to lower pion masses.

Cutoff effects
In general we expect O(a) discretization effects in a HLO µ as the vacuum polarization receives off-shell contributions and the lattice regularization used here is only on-shell improved. Some indications on the size of cutoff effects can be gathered from Fig. 4, where at least for low pion masses discretization effects appear to be rather small and below our statistical errors. However, we believe it is more instructive to look directly at lattice artifacts inΠ (N f ) (q 2 ), as those conceivably depend onq 2 . In Fig. 6 we compare the subtracted vacuum polarization from the ensembles N5 and A3, which have rather similar spatial extensions and pion masses but different values of the lattice spacing. As expected, lattice artifacts mainly distort the function at large values ofq 2 , in the region which contributes little to a HLO µ (see Fig. 3).

Finite size effects
As for the case of cutoff effects we look at finite size effects directly inΠ (N f ) (q 2 ), since different states are expected to contribute to the vacuum polarization whenq 2 is varied (see also [40]). In Fig. 7 (left panel) we compare the results from the E5 ensemble with those from a simulation at the same bare parameters but different L/a (D5 ensemble, not in Table 1), specifically L/a = 24 instead of 32. Finite size effects are clearly visible in this rather extreme setup (m π L 3.4 for D5), but they seem to drop below our present statistical error as L is made larger, as the right panel of Fig. 7 shows. Note, however, that in this second comparison both L and a differ.

CONCLUSIONS
We have presented a calculation of the leading hadronic contribution to the anomalous magnetic moment of the muon on the lattice with two dynamical flavours and a quenched strange quark. We have discussed technical improvements, which led to a better determination of the external momentum dependence of the vacuum polarizationΠ (N f ) (q 2 ). Specifically, together with [41], this paper contains on of the first numerical applications of partially twisted boundary conditions to a quantity containing flavour-singlet currents. The approach follows from the theoretical setup devised in [16]. We restrict the computation to the connected part ofΠ (N f ) (q 2 ) only. As it is clear that presently the accuracy of lattice results is not yet at the level required by phenomenology, the main goals of this paper are a precise assessment of the different sources of uncertainty and an estimate of their size. The main results are listed in Eqs. 3.5 and 3.6, where statistical errors only are included. In the following we estimate the uncertainties due to the modelling of the q 2 -dependence of the vacuum polarization, chiral extrapolations, lattice artifacts and finite volume effects. These uncertainties are quantified by changing one ingredient at a time in the conservative analysis procedure, which we have followed in producing our main results. In detail, on top of the statistical error we identify the following sources: • Fitting procedure. We repeat the entire analysis by using fit ansatz C instead of ansatz A everywhere and obtain which suggests an uncertainty well below our statistical errors.
• Chiral extrapolation. By comparing the results in Eqs. 3.5 and 3.6 to those in Eqs. 3.10 and 3.11, we conclude that at the moment this systematic cannot be resolved with our statistical errors.
• Cutoff and finite size effects. As discussed in Sect. 3.3 these effects appear to be small at the volumes, masses and lattice spacings considered here, but a more comprehensive study is required.
• Uncertainty in the lattice spacing. In [25] the lattice spacing at β = 5.3 is given as a = 0.063 (3)  We will include this systematic by taking half the difference with respect to the results in Eqs. 3.5 and 3.6.
In general we see that at the moment most of these systematic effects, with the exception of the uncertainty on a, are well below our statistical errors. They will become relevant and will have to be more precisely estimated once the latter are reduced. We quote as final results the values determined at β = 5.3, having used the other two lattice spacings to estimate systematic uncertainties. No continuum extrapolation has been performed yet. The present computation of the connected contribution to the anomalous magnetic moment of the muon then gives where we have combined in quadrature statistical and systematic errors. The overall error can definitely be reduced, as it is still statistics-dominated. However, once it is evaluated, the contribution from disconnected diagrams, which is estimated to be around 10% (see [16]), will become the main uncertainty. We therefore consider studying the accuracy that can be reached on the numerical estimates of the disconnected contribution as a priority. A combination of the methods discussed in [42][43][44][45] seems very promising in this respect, and we have also started implementing similar techniques in the context of mesonic three-point functions [46]. Further improvements include considering lighter pion masses and enlarging the set of simulations at β = 5.5.