Study of systematic uncertainties in hadronic vacuum polarization contribution to muon $g-2$ with 2+1 flavor lattice QCD

We study systematic uncertainties in the lattice QCD computation of hadronic vacuum polarization (HVP) contribution to the muon $g-2$. In this article we investigate three systematic effects; finite volume (FV) effect, cutoff effect, and integration scheme dependence. We evaluate the FV effect at the physical pion mass on two different volumes of (5.4 fm$)^4$ and (10.8 fm$)^4$ using the PACS10 configurations at the same cut-off scale. Our results indicate that the chiral perturbation theory (ChPT) possibly underestimates FV effects in a long distance region, where multi-hadron state contributions including two-pion state become prominent. For the cutoff effect, we compare two types of lattice vector operators, which are local and conserved (point-splitting) currents, by varying the cutoff scale on a more than (10 fm$)^4$ lattice at the physical point. The local-local current correlator shows smaller scaling violation than the local-conserved one both in light and strange quark channels. For the integration scheme dependence, we compare the results between the coordinate- and momentum-space integration schemes at the physical point on a (10.8 fm$)^4$ lattice. Our result for the HVP contribution to the muon $g-2$ is given by $a_\mu^{\rm hvp} = 737(9)(^{+13}_{-18})\times 10^{-10}$ in the continuum limit, where the first error is statistical and the second one is systematic. %This value is compared with the phenomenological estimate and recent lattice QCD results.


Introduction
The muon anomalous magnetic moment (g − 2) µ has been a key observable for a proof of predictability of quantum field theory. We expect that there might be a sign of the new physics beyond SM (BSM) in the muon g−2 anomaly, which is 3-4 sigma deviation between the standard model (SM) prediction and the BNL experiment [1,2] suggested in 2004. In order to establish that the (g − 2) µ experiments in FermiLab and J-PARC [3,4] aiming at a factor of 4-5 improvement from the BNL experiment is forthcoming. However, the high precision experiments are not sufficient for the search of the BSM physics [5] since the magnitude of theoretical uncertainty in the SM prediction has not been comparable to that in the new experiments yet. The biggest uncertainty left in the SM prediction is coming from the hadronic vacuum polarization (HVP) effect, which is the leading order of the hadronic contribution to (g − 2) µ denoted by a hvp µ . The phenomenological estimate of a hvp µ [6][7][8][9][10][11], which has been employed in the SM prediction, is obtained by the integrated hadronic R-ratio measured in e + e − annihilation experiment. Including several hadronic decay channels with a particular choice of the center-of-mass energy √ s window, in which the perturbative QCD is used for √ s ≃ 2 GeV, a hvp µ is phenomenologically estimated at 0.4% level of precision [12].
Lattice QCD (LQCD) is another approach to estimate a hvp µ totally independent of the phenomenological estimate. This is a theoretical calculation based on the first principle of QCD, whereas the current precision of LQCD estimate is roughly an order of magnitude larger than the phenomenological one, and it then does not satisfy the accuracy to search the BSM physics (see a recent review [13] and references therein). The main difficulty of LQCD calculation is that, in the Euclidean space-time, the detailed behavior of the HVP contribution with high precision are required around the peak position of the QED kernel, which is much below the hadronic scale of Λ QCD . In such a low-energy region, corresponding to a long distance in the coordinate space-time, it is not an easy task to make a high precision measurement of the HVP contribution because of the exponentially diminishing signal-to-noise ratio in a deeply infrared regime. In addition, the contribution of the ρ resonance state decreases in this regime, while two-pion or three-pion state contributions, which are possible decay modes of the vector resonance, become prominent. This means that a sufficiently large volume at the physical point, where the ρ resonance has open-threshold and the multi-pion states are allowed, is required in the LQCD calculation to correctly estimate the HVP contribution. Furthermore, it is imperative for LQCD to take account of the cutoff effect to obtain a hvp µ in the continuum limit. So the LQCD determination of a hvp µ at a sub-percent precision is still a challenging task. Recent LQCD calculations [14][15][16][17][18] are carried out with the aid of an estimate of effective models, for instance, the chiral perturbation theory (ChPT) [19][20][21] or the Gounaris-Sakurai (GS) parametrization [13,15], to correct the FV effect on ( ∼ < 6 fm) 3 boxes at a long distance; in Ref. [17] the leading-order ChPT estimate is added to the lattice result on a (5.4 fm) 3 box at the physical pion mass taking higher order contributions of O(p 4 ) as a systematic error; Ref. [16] employs a similar strategy to add the ChPT estimate to the lattice results on (6.1-6.6 fm) 3 lattices around the physical pion mass but taking itself as a systematic error conservatively; in Ref. [15] the GS parametrization is used to fit the LQCD result of the vector correlator on a roughly (4 fm) 3 lattice at the unphysical pion mass (m π ≥ 185 MeV) with the time-slice cut of 1.1 < t cut < 1.4 fm; Refs. [14,18] take account of only the two-pion contributions based on the analytic estimate with ChPT.
As pointed out in our previous study [22] it is essentially important to assess the FV effect in the LQCD calculation of a hvp µ employing the direct comparison between different volumes at the physical pion mass without any reliance on the effective models. We have made the direct evaluation of the systematic uncertainty of the FV effect using different volumes of (8.1 fm) 3 and (5.4 fm) 3 near the physical pion mass (m π =135-145 MeV). The difference between the results on two volumes was found to be larger than the ChPT estimate, though the statistical error was so large that they are consistent within 1 σ error bar. In this article, we perform more precise comparison with ChPT using an enlarged lattice more than (10 fm) 4 at the physical pion mass, which are a subset of the PACS10 configurations [23] generated by the PACS Collaboration. We also investigate the lattice cutoff effect by comparing the results at two different cutoffs of a −1 = 2.33 and 3.09 GeV keeping the physical volume at more than (10 fm) 4 . We finally estimate an extrapolated value of a hvp µ in the continuum limit and compare it with other recent LQCD results. This paper is organized as follows. In Sec. 2 we explain the notation and the LQCD methodology to calculate a hvp µ . Lattice parameters and numerical method are explained in Sec. 3 is given by the integral of the vacuum polarization function (VPF) Π(Q 2 ) from zero to infinity in terms of the space-like momentum squared Q 2 M = −Q 2 < 0: whereΠ(Q) is scheme-independent due to a subtraction of scheme-dependent Π(0). The QED kernel K E (s), which is obtained by the one-loop perturbation with α e = 1/137.03599914 and m µ = 105.6583745 MeV [24], has a sharp peak at Q 2 ≈ ( √ 5 − 2)m 2 µ = 0.003 GeV 2 and a rapid falloff for Q 2 → 0. The vacuum polarization function can be extracted from a factorization of the vacuum polarization tensor Π µν (Q), which is given by the Fourier transformation of the vector-vector current correlator, where the index Γ in the superposition of the vector current V µ denotes two choices of lattice operators. One is the local current with Γ ′ = L: with Z V the renormalization constant, and the other is the conserved current with Γ = C: in the point-splitting form with the link variable U µ (x), which preserves the lattice Ward-Takahashi identity µ ∇ * µ V µ = 0 with the backward differential ∇ * (x, y) = δ x,y − δ x−μ,y in the naive Wilson quark action. Note that the lattice local vector currents are not O(a) improved in this study, since we find negligibly small scaling violation (see Section 4.2).
The expression of Π µν (Q) in Eq. (5) has extra contributions of O((aQ) n ) with n ≥ 2 due to the Lorentz symmetry breaking on the discretized space-time in LQCD. After subtracting these lattice artifacts [25][26][27]Π(Q 2 ) computed with LQCD is consistent with the perturbative representation of Adler function [28] in high Q 2 > 1 GeV 2 except for the non-perturbative objects such as d-dimensional operator condensate term given by O d /Q 2d appearing in the operator product expansion (OPE) [29]. For the actual computation of a hvp µ , the LQCD evaluation of the integral of Eq. (1) can be replaced by the perturbative one in the high Q 2 region from some particular point of Q 2 pQCD to infinity. Practically, the integrand for Q 2 pQCD > 1 GeV 2 in Eq. (1) gives minor contribution to total a hvp µ so that the OPE contribution should be negligible. We will discuss it later.
In LQCD we need to evaluate Π(0) by the extrapolation of VPF to the zero momentum limit. Since the minimum momentum in LQCD is defined as Q min = 2π/L, a large volume allows us to perform the qualified extrapolation with less uncertainty of fitting procedures. Once Π(0) is determined, the momentum integral of Eq. (1) is straightforwardly performed with the extrapolation function in the low energy regime, and we can add the perturbative QCD formula in high energy regime of Q 2 > Q pQCD . As pointed out above, the choice of Q 2 pQCD >1 GeV 2 gives only minor contribution to the total a hvp µ . In our analysis the Q 2 integral of Eq. (1) is split into the fit region, the lattice data region, and the pQCD region: where we define lattice momentum Q n = 2πn/L with integer n = {0, 1, ..., L/a − 1} andΠ pQCD is an analytic form in pQCD.Π f (s) ≡ Π f (s) − Π(0) is a functional form with the fitting ansatz, which is used for the extrapolation of the lattice data to obtain Π(0). We utilize three types of fitting functions, (Linear approx.) Π linear f (s) = Π(0) + sY (13) with the fitting parameters Π(0), X 0 , X 1 X 2 . Note that we use the same form of Padé approximation as in Refs. [15,30].

Coordinate-space integration scheme
As an alternative approach we consider the vector-vector current correlator in the coordinate space: where the summation over the same component of sink and source vector currents is taken. With the use of C(x), a hvp µ can be expressed as follows: 1 1 See Appendix A for derivation.
Here there is a degree of freedom for a choice of four components in the momentum Q µ satisfying Q 2 = ω 2 . If we take one component is non-zero and others are zero, i.e. Q µ = {Q ρ = |Q|, Q µ =ρ = 0}, ρ = {x, y, z, t}, Eq. (15) is simplified as where we defineC This is a formula called time-momentum representation (TMR) [31], in particular, with the choice of the time direction for x ρ 2 . It is consistent with the Lorentz-covariant coordinate space representation [32] if the coordinate space integral in Eq. (15) is transformed into the spherical and radial integrals.
In LQCD we perform a discretized coordinate-space summation of the correlator on the finite lattice volume defined as with where r, which denotes a distance from the source point, is regarded as the generalized expression of x ρ in Eq. (16) and thus C ΓΓ ′ (r) representsC ρ (x ρ ) on the lattice. This procedure introduces the systematic uncertainties due to not only the discretized summation but also the truncation at some finite distance r cut 3 . As we will explain below, the lattice used in this study is symmetric and its spatial/temporal extension is large enough to control the finite volume effect and the backward state propagation (BSP) effect investigated in Ref. [22]. We can perform the integral of Eq. (16) (summation of Eq. (18)) for each direction of ρ = x, y, z, t, which allows us to increase the statistics by four times without much computational cost.

Calculation details 3.1 Configurations
We use two subsets of the PACS10 configurations, which are generated with the stout-smeared O(a)-improved Wilson-clover quark action and Iwasaki gauge action [33] on 128 4 and 160 4 lattices (spatial extension L and temporal extension T are symmetric) at β = 1.82 and 2.00, respectively. In addition we also employ the gauge field configurations on a 64 4 lattice at β = 1.82, which are copied in the temporal direction extending T /a to 128 in the FV study. Lattice parameters for these configuration sets are summarized in Table 1. We investigate the FV effect using the 128 4 and 64 4 lattices at the same lattice spacing, and the study of the cutoff effect uses the 128 4 and 160 4 lattices with the fixed physical volume.
The detailed description of the configuration generation on the 128 4 and 64 4 lattices was already given in Ref. [23]. Here we explain the configuration generation on the 160 4 lattice at β = 2.00. We employ the stout smearing parameter ρ = 0.1, and the number of the smearing steps is six, which are the same as in the case of the 128 4 lattice at β = 1.82 [23].The improvement coefficient of c SW = 1.02 is nonperturbatively determined by the Schrödinger functional (SF) scheme following Ref. [34]. The hopping parameters for the light (degenerate up-down) and strange quarks (κ ud ,κ s )=(0.125814,0.124925) are carefully adjusted to yield the physical pion and kaon masses (m π ,m K )=(135.0 MeV,497.6 MeV) with the use of the cutoff of a −1 = 3.09 GeV (a = 0.064 fm) [35] determined from the Ξ mass m Ξ = 1.3148 GeV. The hopping parameter for the charm (only valence quark) is set to κ c = 0.110428 on 128 4 lattice, and κ c = 0.11452 on 160 4 lattice adjusted to physical point.
The degenerate up-down (ud) quarks are simulated with the domain-decomposed HMC (DDHMC) algorithm [36,37] on the 160 4 lattice. The ud quark determinant is separated into the UV and IR parts after the even-odd preconditioning. We also apply the twofold mass preconditioning [38,39] to the IR part by splitting it intoF IR , F ′ IR and F ′′ IR . This decomposition is controlled by two additional hopping parameters: κ ′ ud = ρ 1 κ ud with ρ 1 = 0.9997 and κ ′′ ud = ρ 1 ρ 2 κ ud with ρ 2 = 0.9940.F IR is derived from the action preconditioned with κ ′ ud . The ratio of two preconditioners with κ ′ ud and κ ′′ ud gives F ′ IR . F ′′ IR is from the heaviest preconditioner with κ ′′ ud . In the end the force terms consist of the gauge force F g , the UV force F UV and the three IR forces F ′′ IR , F ′ IR andF IR . The IR forces are obtained with the mixed precision nested BiCGStab method for the quark solver [40]. We adopt the multiple time scale integration scheme [41] in the molecular dynamics (MD) steps. The associated step sizes are controlled by a set of in- Our choice of (N 0 , N 1 , N 2 , N 3 , N 4 ) = (8, 2, 2, 2, 20) for the 160 4 lattices results in 82% acceptance rates. The strange quark is simulated with the RHMC algorithm [42] choosing the force approximation range of [min,max]=[0.000190,1.90] with N RHMC = 10 and δτ s = δτ ′′ IR for the step size. The renormalization constant Z V for the local vector current operator in Eq. (6) depends on the lattice cutoff scale. We obtain Z V = 0.95153(76) at a −1 = 2.33 GeV (β = 1.82) with the SF scheme [43], and Z V = 0.9673 (19) at a −1 = 3.09 GeV (β = 2.00) from the nucleon form factor. Note that we observe a good consistency between the results of Z V determined by the SF scheme and the nucleon form factor [44]. The physical observables are measured at every 10 trajectories on 128 4 and 64 4 , and every 5 trajectories on 160 4 . The statistical error is estimated by the jackknife analysis with 1, 4, and 5 jackknife binsizes for the 128 4 , 160 4 , and 64 4 lattices respectively [23].

AMA with deflated SAP preconditioning
The precision of the light flavor vector-vector current correlator in the infrared (IR) regime, which is the small Q 2 region in Eq. (1) or the long distance region from the source location in Eq. (18), has a vital importance to achieve a less than 1% level of accuracy for a hvp µ with LQCD. As in the previous study [22] we utilize the optimized all-mode-averaging (AMA) technique [45][46][47] to make an efficient calculation of the vector-vector current correlator in LQCD. For the AMA approximation [22,47], we use the parameter set illustrated in Table 2. As shown in Refs. [22,47], the combination of AMA with the deflated Schwartz Alternative Procedure (SAP)  preconditioning [48] achieves the remarkable performance on the large lattice, and it then allows us a precise calculation of a hvp µ , especially, in a long distance region. In fact the condition number in the AMA method with the deflated SAP preconditioning does not have large volume dependence [22] so that the computational cost to solve the light quark propagator does not increase even if the lattice size is enlarged. In addition, the generation of the deflated fields takes only a few percent of total computational cost [22], and hence we do not need to keep the deflation fields in a storage unlike the case of the low-lying mode deflation.
In the left panel of Fig. 1 we show the volume scaling for the relative error of the correlator at the physical pion mass. This is more robust test of volume scaling than the previous study [22], where there might be possible contamination due to the pion mass difference between two volumes. From this plot, one can see that the ratio of the relative error between 64 4 and 128 4 lattices has a consistent behavior with the expected scaling value of 64 3 /128 3 in a long distance region t ∼ > 1.5 fm, which means that the use of a large volume can significantly reduce the statistical error, especially for the IR regime. As illustrated in the right panel of Fig. 1, we also observe the universal behavior for the relative error of the vector-vector current correlator at different cutoff scales on the same physical volume. This feature is also expected from the volume scaling hypothesis for the statistical error.

Multi-hadron state contributions
Using our large lattice ensembles at the physical pion mass, the multi-hadron state contributions, mostly the two-pion state, are correctly involved into the vector-vector current correlator. Figure 2 plots the effective mass of the vector-vector current correlator with (Γ, Γ ′ ) =(L,C) on each ensemble. The ρ meson is allowed to decay into the energetic pions on those ensembles, since the two-pion state energy 2 m 2 π + (2π/L) 2 is much lower than the ρ meson mass m ρ = 770 MeV. We observe that the effective mass goes down below 770 MeV around t ≈ 1 fm and stays above the energy level of 2 m 2 π + (2π/L) 2 in the large t region on each ensemble. This is a clear indication of the existence of lower energy state than ρ meson mass in the region of t ∼ > 1 fm, which is dominated by the multi-hadron state contributions.

Numerical results
With the use of the gauge field configurations explained in Sec. 3, we perform a systematic study of uncertainties stemming from the FV effect, the cutoff effect and the integration scheme dependence in LQCD calculation of a hvp µ . For the FV effect we directly compare the results for the coordinate-space integral of Eq. (18) obtained on the L/a = 128 and L/a = 64 lattices at the same cutoff scale of a −1 = 2.33 GeV. The cutoff effect is investigated by calculating the coordinate-space integral of Eq. (18) on the 128 4 and 160 4 lattices keeping the physical lattice volume constant. We also discuss the operator dependence of the cutoff effect for the local and conserved vector currents. Finally we examine the consistency between the coordinate-and momentum-space integration schemes on the 128 4 ((10.8 fm) 4 ) lattice at the physical pion mass. For the latter we extend T /a to 128 by copying the 64 4 lattice in the temporal direction so that we can eliminate the backward propagation state (BPS) wrapping around temporal direction observed in our previous study [22]. We remark that, although the 64 4 lattice configurations are generated at the same hopping parameter as for the 128 4 lattice, the measured pion mass m π = 139 MeV on the 64 4 lattice is slightly heavier than m π = 135 MeV on the 128 4 lattice due to the FV effect [23]. From the right panel of Fig. 3, one can see that the integrand has a clear tendency that the magnitude of integrand increases when L/a is enlarged from 64 to 128. The left panel of Fig. 4 plots the FV effect defined as  which shows that the magnitude is larger than the leading order ChPT having the same sign with the ChPT prediction [19]. In order to clarify the discrepancy we plot the ratio of the FV effect between LQCD and ChPT at each r cut in Fig. 4. One can observe that the LQCD data tends to become larger than the ChPT prediction from r ≈ 1 fm, and this tendency does not change even if r cut increases, though the statistical error becomes larger.

Finite volume effect
Choosing r cut ≃ 2.6 fm, which is the maximum point of the window method in Ref. [17], the discrepancy of FV effect between LQCD and ChPT in the light quark sector is estimated as on L = 5.4 fm at the physical pion mass. Our result indicates that the actual FV effect is possibly much larger than the ChPT prediction, which could impact on other recent LQCD results using ChPT and similar analytic formula to correct the FV effect for a hvp µ on (4-5 fm 3 ) box [14,[16][17][18]. Our LQCD study suggests that the multi-hadron state contributions not captured in ChPT may become important for FV correction to HVP contribution.
The similar analysis is made for the strange quark sector [a hvp µ ] s . Figure 5 shows little FV effect for [a hvp µ ] s as expected from the fact that the strange quark mass is much heavier than the light one.

Cutoff effect
In Fig. 6 we plot C ΓΓ ′ (r)W r (r) in Eq. (18) at two different cutoff scales of a −1 = 2.33 GeV and 3.09 GeV on the same physical volume over (10 fm) 4 at the physical pion mass. We compare the cutoff effect in two types of the vector-vector current correlators with (Γ, Γ ′ ) =(L,L) and (C,L) for the sink and source vector current operators in Eq. (20). We observe that C LL (r)W r (r) at different cutoff scales well agree with each other, whereas the sizable deviation is found for C CL (r)W r (r) from r ≃0.5 fm. Our LQCD results show the C LL (r) correlator has smaller cutoff effect than the C CL (r) one. In order to make a quantitative measurement of the discrepancy between two types of the correlators, we plot the normalized difference defined as   in Fig. 7. The quantity shows a clear deviation from zero, and its magnitude is reduced for the finer lattice. Around r ≈ 1.5 fm we obtain ∆ r (r ≈ 1.5fm) = 0.089(3) (a −1 = 2.33 GeV) 0.063(1) (a −1 = 3.09 GeV) (24) and their ratio ∆  (1). This suggests that the LQCD result with C CL (r) is affected by the O(a) correction due to a significant cutoff effect on the conserved (point-splitting) current.
In Fig. 8 we plot r cut dependence for [a hvp µ ] lat (r cut ) in the light and strange quark sectors. They asymptotically reach constant values around r cut ∼ > 3.5 fm without large statistical fluctuation. In both the light and strange quark sectors, [a hvp µ ] lat in the (L,L) channel at two cutoff scales agree within 1-σ statistical error, while there is 10-11% cutoff effect to affect [a hvp µ ] lat in the (C,L) channel at a −1 = 2.33 GeV.
We summarize the scaling properties for [a hvp µ ] l lat , [a hvp µ ] s lat and [a hvp µ ] c lat at two cutoff scales and their continuum extrapolations in Fig. 9, where the LQCD results at each cutoff scale are obtained by choosing r cut ≈ 3.5 fm. One can see that the (L,L) channel has rather small cutoff effect, which is not significant in the currently statistical precision, compared to the (C,L) channel in the light and strange quark sectors. In our analysis, a constant fit of (L,L) channel is used for [a hvp µ ] l lat and [a hvp µ ] s lat to take the continuum extrapolation. The systematic error is evaluated by taking the maximum difference between the central value obtained by the constant fit and the linearly extrapolated values in the (L,L) channels with the ansatz of O(a) term including the error of the lattice cutoff itself. The magnitude of the systematic error is comparable with that of the statistical one in the light and strange quark sectors. For the charm quark sector the bottom panel of Fig. 9 shows the large cutoff effect due   to the O(am c ) contribution even in the (L,L) channel. So we take the linearly extrapolated value in the (L,L) channel as the central value in the continuum limit, and its systematic error of the O(a 2 ) contribution is naively estimated as (c 1 a) 2 /c 0 , where c 1 is defined in the fit result of the linear extrapolation c 0 + c 1 a. Further analysis of cutoff effect in charm sector will be done by adding the data on one more fine lattice in the future. One can see that the uncertainty of the cutoff effect is dominant in the charm quark sector. For the total contribution of [a hvp µ ] l lat + [a hvp µ ] s lat + [a hvp µ ] c lat , the uncertainty in the light quark sector is still dominant.

Analysis of momentum-space integration scheme on 128 4 lattice
Compared to the coordinate-space integration scheme, the momentum-space integration scheme is rather straightforward to perform the integral in Eq. (1) once Π(0) is determined by the zero momentum extrapolation of VPF (see Eq. (8)). The advantage in this work over prior ones [15,30,[49][50][51] is that we can not only access VPF in the low-momentum region but also have high resolution in terms of Q 2 without resort to the twisted boundary condition for the valence quark by using a large lattice size more than (10 fm) 4 . Our large lattice is also useful to reduce the uncertainty due to zero momentum extrapolation and remove the partially quenching effect due to the different boundary conditions between sea and valence quarks. We plot the LQCD data of VPF in each quark sector in Fig. 10, where the fit results of Padé approximation of order [1,1] for the light and strange quark sectors using Q 2 fit ≤ 0.05 GeV 2 and linear function for charm quark using Q 2 fit ≤ 0.065 GeV 2 are also shown. We observe that VPF in the light and strange quark sectors have stronger slope around the zero momentum region than that in the charm quark sector. As shown in Fig. 10, the Padé approximation of lowest order of [1,1] well describes such a lattice data with reasonable χ 2 /dof< 1 in the correlated fit. Although we employ a narrow fit range close to zero-momentum, Fig. 11 shows that the fit function agrees with the Q 2 dependence of LQCD data up to Q 2 =0.4 GeV 2 , which is far beyond the fitting range. This behavior indicates the lowest Padé approximation, which consists of single pole dominance, is a reasonable approximation in the IR regime. Since VPF multiplied by the weight function W q in Q 2 ≥ 0.5 GeV 2 gives tiny contribution to the total a hvp µ as mentioned in Sec. 2.1, we evaluate the integral without pQCD part (the third term of Eq. (8)) in our analysis. In fact, the LQCD data of integral more than Q 2 ≈ 0.5 GeV 2 is below 0.5×10 −10 corresponding to less than 0.1% for a hvp µ (see the right panel of Fig. 11), and it is then negligible. Therefore we hereafter estimate a hvp µ using integral up to Q 2 = 0.5 GeV 2 . Since the integrand has the sharp peak structure much below the minimum momentum squared Q 2 min ≈ 0.013 GeV 2 allowed on our ensemble (see Fig. 11), the integral in Eq. (8) is sensitive to the extrapolation procedure from Q 2 min to zero. We employ the linear extrapolation and the Padé approximation of order [1,1] and [2,1]. Figure 12 compares the results of [a hvp µ ] Mom obtained by both extrapolation methods with varying the fitting ranges from Q 2 min to Q 2 fit . In case of the linear extrapolation the results of [a hvp µ ] Mom show significant Q 2 fit dependence, due to the higher order term than O(Q 2 ) even in Q 2 ≈ 0.013 GeV 2 , except for the charm quark sector. On the other hand, we observe little Q 2 fit dependence for [a hvp µ ] Mom with the Padé approximation of order [1,1] and [2,1] up to Q 2 fit = 0.235 GeV 2 in the light quark sector. The results with different Q 2 fit and different orders of Padé approximation are in good agreement with [a hvp µ ] CL lat (r cut ) with r cut = 3.5 fm within 1.5 σ error. We notice the strong Q 2 fit dependence for the results in the strange quark sector appears in Fig. 12. In this case, an extra lattice cutoff effect of O(am s (aQ) 2 ), which is not described by the naive Padé approximation, may arise in strange quark sector. More detailed study will be needed in the future.
In contrast to the 128 4 lattice, Fig. 13 shows the significant Q 2 fit dependence for the results with both extrapolation methods on the 64 4 lattice since the low Q 2 data has coarse resolution on this lattice. Our LQCD study suggests that the lattice size with L = 5.4 fm at the physical pion mass, corresponding to m π L = 3.8, is not large enough for the momentum-space integration scheme to obtain a reliable result of a hvp µ because of large FV correction. We remark that the statistical precision of the result for [a hvp µ ] Mom is more easily obtained than that for [a hvp µ ] lat . This is because of a noise cancellation inΠ(Q) between the extrapolated Π(0) and Π(Q), which are highly correlated with each other. In addition, in the momentumspace integration scheme, we do not need to introduce the truncation of the integration range corresponding to the IR truncation r cut in the coordinate-space integration scheme. It indicates a possibility that once we have low Q 2 data covering a peak position of W q (Q 2 ∼ 0.003 GeV 2 ), for which we need to prepare twice larger box size than this study, we can obtain a high precision result with less statistical and systematic errors than the coordinate-space integration scheme.

Discussion
We obtain the connected a hvp where the first error is statistical for (Γ, Γ ′ ) =(L,L) with the constant fit and the second one is systematic for the uncertainty in the continuum extrapolation explained in Sec. 4.2. We find that the statistical and systematic errors for the light quark sector gives the leading contribution to the total error. The contributions from the strange and charm quark sectors are minor effects.
Here we have two remarks; 1. Our choice of r cut ≈ 3.5 fm in the coordinate-space integration scheme, which is larger than 3 fm employed in Refs. [16,17], is large enough to control the IR truncation. In Figs. 6 and 8 we observe that the integrand has a nonzero value of 23(10) × 10 −10 at r cut ≈ 3 fm in the (Γ, Γ ′ ) =(L,L) channel on the 128 4 lattice and the integral is still increasing, while the integrand is consistent with zero at r cut ≈ 3.5 fm and the integral does not depend on r cut even if we use larger r cut . High precision data on the large lattice more than (10 fm    at the physical point allow us to evaluate the integral with the IR truncation effect under control.
2. The scaling properties presented in Sec. 4.2 are similar to the domain-wall fermion case [17], though the computational cost is much lower for the Wilson-type quark action. The continuum extrapolation is straightforward and theoretically robust for the Wilson-type quark action compared to the staggered fermion case [16,53].
In this paper, we concentrate on the connected HVP diagram, while there are some missing diagrams of the isoscalar contribution with the disconnected diagram and the isospin breaking (IB) term due to the QED correction. Referring to the recent works in Refs. [16,17], we conservatively add the systematic error of the quark disconnected diagram contribution as −2% effect and the IB effect as +1% error to the total contribution. We then find where the first error is statistical and the second one represents the total systematic error obtained in the quadrature. The magnitude of the error is still 2.7%, in which the systematic error, mainly due to the uncertainty of disconnected diagram, is more than twice larger than the statistical one. Compared to other lattice results (N f ≥ 3) (see Fig. 14), our value is consistent with the results by RBC-UKQCD [17] and BMW [16] collaborations, while we find a slight tension from the recent results by ETMC [18], HPQCD [14] collaborations, and two-σ deviation from the phenomenological estimates [10,11]. Our result seems to favor the "experimental" a hvp µ , which is defined by the difference between the BNL experimental value of a µ and the theoretical calculation with QED and EW including the light-by-light scattering contribution in Ref. [12].

Summary
We have studied the systematic uncertainties in the LQCD calculation of a hvp µ on the PACS10 gauge configurations which have more than (10 fm) 4 box size at the physical point with two different lattice cutoffs. This study and previous work [22] are the direct LQCD calculation without use of any ansatz and reliance on any effective models. The optimized LQCD calculation of HVP on sufficiently large lattice size at the physical point allows us to access the deep IR regime where the contributions of multi-hadron states become manifest. Our study reveals that such contributions can not be completely captured by ChPT. In Fig. 14 we observe that our result of a hvp µ is relatively larger than other LQCD studies. The reason of such a tendency may be due to a discrepancy between LQCD and ChPT (or related phenomenological models) including only two-pion state contribution, which were applied to evaluate the FV correction in other LQCD studies, as discussed in Sec. 4.1. We have also investigated the lattice cutoff effect in the coordinate-space integration scheme using the data at two different cutoffs. We find that the cutoff effect is reasonably small for the local vector current on our gauge configurations. Furthermore, the momentum-space integration scheme on L > 10 fm lattice yields the high quality data of VPF close to Q 2 = 0, which substantially reduce the uncertainty in the zeromomentum extrapolation. With a careful study of the extrapolation procedure dependence we have confirmed the consistency between the results in the momentum-and coordinate-space integration schemes.
The total error for the result of a hvp µ is 2.7%, in which the statistical error is 1.2% and the remaining is the systematic uncertainty. We plan to reduce the both statistical and systematic errors with additional calculation including one more finer lattice, disconnected diagram, and   [16], ETMC [18], HPQCD [14], RBC-UKQCD [17] collaborations and phenomenological estimate obtained with the experimental R-ratio by DHMZ [10] and KNT [11]. Shaded vertical band shows the "experimental" a hvp µ estimated by the difference between the BNL experimental value of a µ and the theoretical value with QED and EW including the light-by-light scattering contribution. The error bar for [a hvp µ ] l in this work represents the combined error with the statistical one and the systematic one due to cut-off effect. Additional uncertainties of missing disconnected diagram and IB effect are included in the error bar of a hvp µ in this work.
QED effect in future. Here we will point out a possibility that the momentum-space integration scheme with L > 20 fm covers the peak position of kernel function in low Q 2 regime so that it could be a rigorous test for the LQCD scheme. We leave it to future work.

A The derivation of coordinate representation
Fourier transformation of Eq. (14) is defined as Π in Eq. (2) can be represented aŝ where the second term is expanded as and we obtainΠ In general, the second derivative of G with respect to ω = |Q| can be expressed as where the terms with odd power of x µ vanish in the coordinate integral. Substituting the above equation into Eq. (1), we can obtain Eq. (15).