Hadronic vacuum polarization in the muon $g-2$: The short-distance contribution from lattice QCD

We present results for the short-distance window observable of the hadronic vacuum polarization contribution to the muon $g-2$, computed via the time-momentum representation (TMR) in lattice QCD. A key novelty of our calculation is the reduction of discretization effects by a suitable subtraction applied to the TMR kernel function, which cancels the leading $x_0^4$-behaviour at short distances. To compensate for the subtraction, one must substitute a term that can be reliably computed in perturbative QCD. We apply this strategy to our data for the vector current collected on ensembles generated with $2+1$ flavours of O($a$)-improved Wilson quarks at six values of the lattice spacing and pion masses in the range $130-420\,$MeV. Our estimate at the physical point contains a full error budget and reads $(a_\mu^{\rm hvp})^{\rm SD}=68.85(14)_{\rm stat}\,(42)_{\rm syst}\cdot10^{-10}$, which corresponds to a relative precision of 0.7\%. We discuss the implications of our result for the observed tensions between lattice and data-driven evaluations of the hadronic vacuum polarization.

The anomalous magnetic moment of the muon, as a low-energy precision observable, has long served as a test of the Standard Model (SM) of particle physics.Its experimental measurement has reached a spectacular precision.With the latest result of the Fermilab (g−2) experiment at the 0.20 ppm level [1], the stakes for the SM-based theoretical prediction have been raised even further.The 2020 White Paper (WP2020) of the Muon (g −2) Theory Initiative [2] arrived at a result with a quoted precision of 0.37 ppm, the uncertainty being dominated by the hadronic vacuum polarization (HVP) contribution and to a lesser extent by the hadronic light-by-light contribution.The very recent experimental result [1] is in tension with the 2020 prediction at the level of five standard deviations.However, several results and observations of the last three years concerning the HVP contribution have cast doubt on the reliability of the 2020 prediction [2].
First, while the WP2020 evaluation of the HVP contribution a hvp µ was based solely on the dispersion relation relating it to the inclusive e + e − → hadrons cross section, a lattice QCD calculation with subpercent precision published by the BMW collaboration in 2021 [3] arrived at a result larger by 2.1 combined standard deviations than the WP2020 one.If one were to replace the WP2020 evaluation of the HVP contribution by that of the BMW collaboration, the tension between the SM and the latest experimental result for the muon (g − 2) would shrink to 1.7 σ.
Thirdly, a new measurement of the e + e − → π + π − cross section by the CMD-3 collaboration has appeared, which lies higher than the previous measurements.This new result is very significant, since in the dispersive evaluation of a hvp µ , the π + π − channel dominates, with the ρ meson region 600 MeV < √ s < 900 MeV providing more than half of a hvp µ .If confirmed, the CMD-3 result would bridge the difference between the WP2020 and the BMW calculations of a hvp µ .The 'intermediate window' alluded to above is the second of three subcontributions in a partition of a hvp µ based on the Euclidean time separation between two electromagnetic currents [13,14].The first of these subcontributions is the short-distance quantity (a hvp µ ) SD extending up to 0.4 fm (with a soft edge), which is the focus of this paper, while the third, numerically dominant part is the long-distance contribution (a hvp µ ) LD .The short-distance quantity only represents about ten percent of a hvp µ ; in this respect, its determination at a precision level of two percent would currently be sufficient.However, we argue below that it is worth aiming at a subpercent precision for (a hvp µ ) SD , in order to provide a clue as to the origin of the well established discrepancy in the intermediate window quantity.
If one assumes that the tension between the lattice and the dispersive determinations of the intermediate window is due to an underestimate (by an overall factor of about 0.94) of the experimentally determined R(s) ratio in the interval from 600 to 900 MeV, then one expects to find an underestimate of the dispersively determined short-distance quantity by more than one percent.Specifically, based on our own previous calculation [4] and the estimates given therein of the fractional contributions to (a hvp µ ) SD and (a hvp µ ) ID associated with the interval 600 MeV < √ s < 900 MeV, the expectation would be an excess of (1.44 ± 0.37)% of our lattice result for (a hvp µ ) SD over the dispersive one.Thus it is clear that a precision well under one percent is required on (a hvp µ ) SD in order for such a calculation to have an impact.The existing lattice results [5,6] for (a hvp µ ) SD are consistent with the scenario above; at the same time, the evidence for a discrepancy with the dispersive value for (a hvp µ ) SD lies only at the 1.4 σ level even for the most precise result [5].
There are further hints in favour of the scenario of an underestimated R(s) ratio in the ρ region.One is the calculation of the hadronic contribution to the running of the electromagnetic coupling α(q 2 ).The Mainz/CLS publication of 2022 [15] showed that a tension develops between photon virtualities of zero and 1 GeV 2 , while the running is consistent with the dispersive approach at higher q 2 ; see [3,16] for the earlier calculation by the BMW collaboration.This observation points to an origin at √ s ≲ 1 GeV in the dispersive approach, or alternatively to an issue with the lattice data at Euclidean times x 0 ≳ 1 fm.
Secondly, the fact that the absolute difference in (a hvp µ ) ID appears to reside entirely in the light-quark connected contribution [12] can be interpreted as a further consistency check, since the on-shell effects of the strange quark start at √ s = 2m K ≃ 990 MeV.Thirdly, a dedicated study of the spectral function, smeared by a Gaussian kernel, has recently been published [17], where the strongest tension (about 2.9 σ) between the lattice and the e + e − data-based result is seen around the ρ peak with a broad smearing kernel of width 630 MeV.Last but not least, an increase of the phenomenological R(s) ratio by six percent in the interval 600 MeV < √ s < 900 MeV would bring the WP2020 prediction for the muon (g − 2) into perfect agreement with the direct experimental measurement.For all these reasons, we consider such a scenario to be a useful working hypothesis.At the same time, in this scenario a very convincing explanation (within or outside the Standard Model) would have to be found why all previous experimental measurements of the π + π − channel and possibly also π + π − π 0 channel are systematically suppressed in the ρ, ω region, well beyond what the quoted uncertainties would allow for.The short-distance observable (a hvp µ ) SD can be obtained very precisely in lattice QCD as far as statistical errors are concerned.Some sources of uncertainty that are very relevant to the calculation of a hvp µ , such as finite-volume effects or a slight mistuning of the light quark masses, play hardly any role in (a hvp µ ) SD .Also the sensitivity to the 'absolute scale setting', meaning the calibration of the lattice spacing in physical units, is very subdominant.The main source of systematic uncertainty is the presence of enhanced cutoff effects, which must be removed via a careful continuum extrapolation if one targets subpercent precision.Given that cutoff effects are enhanced at short distances, we devise a specific observable that we subtract from (a hvp µ ) SD to facilitate the continuum limit, and evaluate this observable with the help of massless perturbation theory to order α 4 s .It is designed to only receive contributions from virtualities of order Q 2 , where Q 2 is an adjustable parameter chosen to lie in the perturbative regime.
This paper is organized as follows.In section II, we describe our notation and computational strategy, our gauge ensembles and fitting procedures.Section III contains our results for the various contributions computed using lattice QCD, with two further, small effects evaluated in its last subsection III H.We present our final result in section IV, and discuss its impact on the current (g − 2) puzzle.

II. PRELIMINARIES
Here we recall some basic definitions and describe our overall computational strategy.In particular, we introduce and discuss the decomposition of the short-distance window observable, which allows us to gain good control over the extrapolation to the continuum limit through a combination of perturbation theory and lattice calculations.
The electromagnetic current correlator in the time-momentum representation (TMR) is defined by where J µ (x) =2 3 ūγ µ u−1 3 dγ µ d−. . .In this representation, the hadronic vacuum polarization contribution to the muon anomalous magnetic moment is given by 1 The function K(z) is given analytically in Ref. [18].Its asymptotic properties are The short-distance window contribution is given by where is a smooth step function interpolating around x 0 ≈ d between zero and unity on a distance scale ∆, and we choose the currently standard values For the perturbative evaluation of high-momentum contributions to a hvp µ , it is useful to introduce the Adler function, defined as the logarithmic derivative of the vacuum polarization function, Note that D(Q 2 ) is normalized so as to have the same high-energy asymptotics as the R ratio, will play an important role in the following.It can be computed in the Euclidean theory via the integral [19] Π Note that the TMR kernel to compute this quantity is positive-definite, proportional to x4 0 at short distances and remains bounded at long distances.
A. Flavour structure of the current-current correlator We write the current with the help of a matrix T a acting in flavour space, and introduce the flavour-specific TMR correlator For a = 1, . . .8, we define T a = λ a 2 ⊕ 0 c ⊕ 0 b to be given by the corresponding Gell-Mann matrix λ a acting in the (u, d, s) sector.We also use J c µ = cγ µ c (J b µ = bγ µ b) to denote the charm (bottom) current, i.e.T c = diag(0, 0, 0, 1, 0) and T b = diag(0, 0, 0, 0, 1).The index γ is used to indicate the physical quark charge matrix, ).With this flexible notation in place, we can write the electromagnetic-current correlator as As the notation indicates, for the heavy quarks we treat separately quark-connected and disconnected contributions.The dots stand for contributions that are too small to be of relevance in this work, namely disconnected diagrams involving the bottom quark and the top contribution.We use the corresponding notation for the different flavour contributions to a hvp µ ≡ a γ,γ µ as well as for (a hvp µ ) SD ≡ (a γ,γ µ ) SD , in particular

B. General computational strategy
For a contribution to (a hvp µ ) SD of a given flavour, we introduce the decomposition2 where On the basis of eq. ( 3), the (leading) x 4 0 part of the weight function defining (a hvp µ ) SD is cancelled in the integral of eq.(17).We expect this cancellation to reduce cutoff effects in the lattice calculation.On the other hand, the subtracted quantity can be computed using the Adler function exclusively at virtualities greater or equal to Q 2 /4.The known good convergence properties of D(Q 2 ) for Q ≳ 2.5 GeV lead us to expect similarly good properties for b (d,e) .In eq. ( 16) applied to the isovector channel, we intend to compute b (d,e)  by inserting the massless perturbative Adler function into eq.( 9).In other channels, the finite quark-mass effects in b (d,e) can be computed in lattice QCD up to the charm quark mass.
On the left hand side of Figure 1 the kernel function K SD sub (m µ , Q, x 0 ) of eq. ( 18) is shown for several choices of the virtuality Q.The very-short distance region of the integrand for (a hvp µ ) SD is excluded in all cases.In the following, more specific aspects of our strategy are discussed for the most important channels.

The isovector contribution: test of massless perturbation theory
We can write This equation can serve as a test to compare lattice results for the left-hand side to the perturbative Adler function based calculation of the right-hand side.For orientation, at leading order in perturbation theory b (3,3) 2. Strategy for the (u, d, s) isoscalar contribution Knowing (a 3,3 µ ) SD from the lattice, one method of obtaining (a 8,8 µ ) SD is via the decomposition (a 8,8 µ ) SD = (a 3,3 µ ) SD + ∆ ls (a µ ) SD .
Only the term ∆ ls (a µ ) SD needs to be computed anew.It is SU(3) f breaking and therefore parametrically suppressed at short distances: in the continuum, the linear combination of correlators G (8,8) −G (3,3) has the leading parametric behaviour α s (m 2 s −m 2 l )/|x 0 | at |x 0 | → 0, thus making the corresponding integrand for (a hvp µ ) SD very suppressed at short distances.Therefore no help from perturbation theory is needed to compute ∆ ls (a µ ) SD .
Note that the estimator (21) is entirely equivalent to proceeding via eq.( 15), with b (8,8) (Q 2 ) evaluated via b (8,8) = b (3,3) + ∆ ls b , where b (3,3) is taken from massless perturbation theory and ∆ ls b from the lattice.However, it appears simpler to focus on ∆ ls (a µ ) SD , since the latter quantity is expected to be much smaller than the already computed (a 3,3 µ ) SD .

Strategy for the charm-connected contribution
Similarly, the charm-connected contribution can be estimated according to eq. ( 15), where the subtraction function b where the first term is taken from massless perturbation theory and ∆ lc b is evaluated on the lattice.

C. Lattice setup
We refer to our previous publications [4,20] for an in-depth overview of our computational setup.Here, we only point out the refinements that we have implemented in view of computing the short-distance contribution and a hvp µ .

Gauge ensembles
We perform our computation on the 2+1 flavour CLS ensembles with O(a) improved Wilson fermions and a tree-level O(a 2 ) improved Lüscher-Weisz gauge action [23,24].Twistedmass and RHMC determinant reweighting are used to stabilize the simulations and to simulate the strange quark, respectively [25][26][27][28].Compared to our recent computation of the intermediate-distance contribution in [4], we have extended the set of gauge ensembles that is used in our study.By including a second ensemble with physical quark masses we are able to more tightly constrain mass-dependent cutoff effects which otherwise could spoil our strategy to control the continuum extrapolation with ensembles at larger-than-physical light quark masses.
As the sum of the bare sea quark masses is held constant for each value of the bare coupling on the ensembles that have been used so far, the combination of meson masses π is approximately constant along each chiral trajectory.To correct for a small deviation from m phys K when approaching m phys π , due to cutoff effects and higher orders in chiral perturbation theory, we have performed dedicated measurements to compute the derivatives of our observables with respect to the quark masses in [4].This allowed us to constrain the dependence of (a hvp µ ) ID on m 2 K + 1 2 m 2 π .In view of the large statistical uncertainties of these derivatives when considering the region of large Euclidean times, we have adapted our strategy by including four ensembles on a different chiral trajectory where the strange quark mass is approximately held at its physical value.Since the pion masses on these additional ensembles are in the range between 200 MeV and 260 MeV, the deviation from (m phys K ) 2 + 1 2 (m phys π ) 2 is at the level of a few percent on these ensembles.An overview of the ensembles used in this work is given in table I.  [21,22], the pion and kaon masses, the physical size of the lattice and the length of the Monte Carlo chain in Molecular Dynamics Units (MDU).Ensembles with an asterisk are not included in the final analysis but used to control finite-size effects.Ensembles with a dagger are only on the chiral trajectory where m s ≈ m phys s .Ensemble N452, marked with a plus, is only used in the computation of isospin breaking effects.

Id
β bc

O(a) improvement
We perform full O(a) improvement of the action [29] and the observables used in this work.We use the local (L) and the point-split (C) currents on the lattice, defined via where U µ (x) is the gauge link in the direction μ associated with site x.With the local tensor current defined as with two independent sets of non-perturbatively determined coefficients c (α) V from [30] and [31].In past work, we have employed the symmetric discrete derivative ∂ν ) to compute the time derivative of the tensor current.At very short Euclidean distances, where the vector-tensor correlator falls off steeply proportional to x 2 0 , the cutoff effects from the discrete derivative can be significant, despite being O(a 2 ).By rewriting the derivative acts on a function that varies more slowly and the size of the cutoff effects at short distances is reduced.Note that a similar effect may be achieved by rewriting the TMR integral via integration by parts such that the time derivative only acts on the functions K(m µ x 0 ) and w SD (x 0 ) which are defined in the continuum.In this work, we use the derivative of eq. ( 26) to reduce the cutoff effects in the short-distance region.

Tree-level improvement
As pointed out in section II B, we reduce cutoff effects from very short Euclidean distances by evaluating a part of (a hvp µ ) SD perturbatively, thus cancelling potentially dangerous effects of O(a 2 log(a)).A further reduction of cutoff effects from short distances may be achieved by computing these to tree-level in perturbation theory [5,19] and subtracting them from the non-perturbatively calculated observable.Denoting the tree-level evaluation of an observable O(a) by O tl (a), we can improve the approach to the continuum limit by either of the replacements For (a hvp µ ) SD sub , we find that the multiplicative improvement seems to reduce the cutoff effects slightly better than the additive version and therefore use it in this work.Note that the additive improvement may be used to modify a hvp µ via We use massless perturbation theory at leading order to compute (a hvp µ ) SD,tl sub .To all orders in perturbation theory, the ⟨Σ µν (x)J λ (0)⟩ correlator is a pure O(a) artefact in the massless theory.Therefore, the currents are automatically O(a) improved, whereas the tree-level values c L,tl V = 0 and c C,tl V = 0.5 would have to be used for massive correlators.We find that using the non-perturbatively determined values for c (α) V , the same ones as in the nonperturbative computation, leads to a further reduction of cutoff effects.We interpret this as a cancellation of effects of O(a n ), n ≥ 2 that are introduced by the derivative of the tensor current when performing O(a) improvement.

Chiral-continuum fit
Cutoff effects are clearly the biggest challenge for a precise computation of (a hvp µ ) SD from lattice QCD.The reliable extrapolation of our results to the continuum limit with a conservative estimate of the corresponding systematic uncertainty is therefore the main challenge of this work.Even after applying the techniques described above to tame discretization effects, higher order cutoff effects, compared to a 2 scaling, are visible in our data and thus have to be included in our fits.We note that it is not the relative size of the cutoff effects but this curvature that makes the extrapolation challenging.Given the very precise data and the spread of quark masses on the ensembles which are included in our computation, we can also resolve mass-dependent cutoff effects.
We therefore test a variety of ansätze to describe the dependence of our data on the lattice spacing and the quark masses.The fit to the quark mass dependence is no particular challenge since we include two ensembles at physical quark masses and therefore do not need to extrapolate.We find mild quark mass effects in the short-distance region.
We set the scale using the gradient flow scale t 0 /a 2 [32] together with its physical value t phys 0 = 0.1443 (7) fm which was determined in [21] from a combination of f π and f K .As an alternative, we have in the past directly set the scale with the pion decay constant by the ratio f phys π /(af π ).This was found to reduce the slope of both chiral and continuum extrapolations for a hvp µ in [20].Since the short-distance contribution has only a mild dependence on m π that is significantly different from that of f π , we do not use f π -rescaling in this work.However, we have verified that we would obtain very similar results to the ones quoted in this work, albeit with larger uncertainties due to the enhanced difficulty in the extrapolation of the data to the physical point.
We also introduce the dimensionless combinations To describe the chiral dependence of the observables O computed in this work, we always include a term that is linear in m 2 π and allow for another term that includes a higher order correction, where The strange quark mass dependence is parametrized via the parameter which is close to its physical value on all of the ensembles considered in this work.We fit to To describe the cutoff effects with sufficient quality, we need to allow for a number of terms, guided by Symanzik effective theory.Our most general ansatz is where X a = a 2 /(8t 0 ).It is not possible to constrain all fit parameters at once.We build a variety of different descriptions of the chiral and cutoff effects by setting one or several of the parameters γ i , β i , δ i , ϵ 2 to zero.We note that, in an O(a) improved theory, the prediction from Symanzik effective theory for the scaling is a 2 [α s (1/a)] Γ with Γ the leading anomalous dimension [33,34].Any curvature in a 2 may therefore be due to the modified scaling behaviour and not due to higher order cutoff effects. 3The leading anomalous dimension for local quark bilinears with vector quantum numbers has been computed to be Γ = 0 in the O(a) improved action used in this work [35].The lowest anomalous dimension from the gauge action is 0.76 [34] and larger than the smallest non-zero contribution from the current, which is 0.395.We allow for effects from beyond leading order terms by including the modification in our continuum extrapolations, thereby doubling the number of fit ansätze.We use fiveloop running [36], starting at Λ MS from [37], to determine the running-coupling constant.Generally, we find that fits with Γ = 0 are favoured over fits with Γ = 0.395.However, the latter have non-vanishing weights in our averages and prefer continuum extrapolated results that are shifted to slightly smaller values.The effect of including a non-vanishing anomalous dimension is found to be significantly smaller than the overall uncertainty in all cases considered in this work.
We extend our explorations of the parameter space for the chiral and continuum extrapolations by performing cuts in the data and repeating the analysis with reduced data sets.We perform cuts in the lattice spacing by removing the coarsest or the two coarsest lattice spacings from our analysis.Furthermore, we perform cuts by excluding all ensembles with m π < 400 MeV, or with m π < 300 MeV respectively in the case of contributions that vanish at the SU(3) symmetric point, from our analysis.  (3,3− G (8,8) ) FIG. 1: Left: Illustration of the modified short-distance kernels according to eq. ( 17) with varying virtuality.Right: Integrands for the computation of (a hvp µ ) SD on ensemble E250 with near physical quark masses at a ≈ 0.064 fm.Smaller contributions have been scaled as indicated in the plot.
To assign and compare fit qualities in view of the change of the number of fit parameters and data points, we apply the Akaike information criterion [38] and the model averaging method from [39].We assign a weight to each fit, where k i is the number of fit parameters and n i is the number of data points in the fit with minimized χ 2 i .The normalization N is such that i w i = 1. 4e obtain the central value and statistical uncertainty of an observable O by a weighted average over all analyses Our estimate of the systematic error associated with the extrapolation to the physical point is given by (δO Statistical uncertainties are determined and propagated using the Γ-method in the implementation of the pyerrors package [41][42][43].

III. RESULTS
To illustrate the behaviour and the relative size of the contributions to (a hvp µ ) SD , we show in the right panel of Figure 1 integrands of the various contributions according to eq. ( 14), i.e., without subtraction, where the smaller contributions are scaled as indicated in the figure to show them on the same scale.The LC discretization is displayed for the quark-connected contributions and the CC discretization for the two quark-disconnected contributions.The charm-connected correlation function contributes significantly to the final result for (a hvp µ ) SD .Except for the two charm-disconnected contributions, the relative statistical uncertainties are very small.
As outlined above, we non-perturbatively compute the subtracted observables (a d,e µ ) SD sub (Q 2 ) according to eq. ( 17).After the combination with b (d,e) (Q 2 ), eq. ( 16), the dependence on the virtuality Q has to cancel.We perform our calculation for several values of Q to explicitly test this.For our final result, we choose Q = 5 GeV as a good compromise between a high enough scale for the perturbative evaluation of the Adler function and a smoothly varying kernel function K SD sub (m µ , Q, x 0 ) on the lattice.If not stated explicitly, all results for the quantities b (d,e) , (a d,e µ ) SD sub and (a hvp µ ) SD are from here on expressed in units of 10 −10 .
A. Perturbative evaluation of b(Q 2 ) in the isovector channel We have used the perturbative coefficients of the Adler function as given in [44] up to α 3 s and taking the O(α 4 s ) term from [45].The effect is to multiply the parton-level prediction by the factor (1 + (α s /π) + 4 n=2 (α s /π) n d n ).We use the perturbative Adler function in the N f = 3 massless theory using Λ MS = 0.338 (12) GeV [46].The QCD beta function was taken into account up to three-loop order included, and we used the corresponding large-Q 2 asymptotic solution for α s throughout, after having checked that using instead the numerical solution for α s makes an entirely negligible difference.As expected, the convergence of perturbation theory is excellent.For instance, at Q = 5 GeV, the set of values of the subtraction term, from the parton-level prediction to the highest-order prediction read b (3,3) (25 GeV 2 ) = {7.422,7.971, 8.037, 8.057, 8.068}. ( Table II lists the perturbative results for several values of Q 2 .Secondly, the uncertainty of Λ MS induces an absolute uncertainty of 0.011 on b (3,3) (25 GeV 2 ).Thirdly, the sea-quark effect due to the strange quark is small.In the estimates above, the strange-quark mass is treated as zero.Evaluating the expression in the N f = 2 massless theory (corresponding to m s = ∞) with Λ MS = 0.330 GeV [46], we obtain The comparison of Eqs. ( 39) and (40) shows that the effect of changing m s = 0 to m s ≃ 100 MeV must be negligibly small, since it is parametrically suppressed by m 2 s /Q 2 .Finally, we have checked that the phenomenologically estimated gluon and quark 'condensates' [47] make a negligible contribution.
Anticipating that the quantity b (3,3) (25 GeV 2 ) is the only perturbative prediction that enters our final result for (a hvp µ ) SD , we take as its full uncertainty the linear sum of the error from Λ MS and the size of the last (O(α 4 s )) perturbative term, b (3,3) (25 GeV 2 ) = 8.068 ± 0.022.( 41)  The perturbatively computed subtraction terms b (3,3) (Q 2 ) and b (c,c) (Q 2 ) (in units of 10 −10 ) for the isovector and the charm correlator respectively.The former is computed using the N f = 3 massless perturbative Adler function to O(α 4 s ) included, the latter with the massive perturbation theory to O(α 2 s ) included, in the expansion with respect to the ratio (m 2 c /Q 2 ) to third order included.The first number in brackets indicates the absolute size of the last O(α n s ) term included in the estimate and gives an indication of the uncertainty related to the convergence of the perturbative series.For b (3,3) , the uncertainty induced by that on Λ MS is indicated as well.µ ) SD sub (25 GeV 2 ).Left: Best fit, according to its model weight, to the data for set 1, LL.The data is corrected for deviation from physical Φ 4 .The coloured lines show the evaluation of the fit function at finite lattice spacing.The black line, together with the grey band, show the dependence on Φ 2 in the continuum limit.Right: Approaches to the continuum limit for four sets of data based on the improvement schemes of set 1 and 2 and the LL and LC discretizations of the current based on a scan over fit models for (a 3,3 µ ) SD sub (25 GeV 2 ).Each line shows the result from one single fit and the opacity of the lines corresponds to the weight of the fit in the model average.Dashed vertical lines indicate the lattice spacings used in this work.
B. Non-perturbative evaluation of (a hvp µ ) SD sub (Q 2 ) in the isovector channel In the short-distance regime, no signal-to noise problem hinders the computation in the isovector channel and the precise lattice data can be integrated by summation.We evaluate (a hvp µ ) SD sub at several values of Q 2 to be able to monitor the convergence of the sum in eq. ( 15) and to investigate whether the choice of Q 2 has an influence on the systematic uncertainties of the lattice evaluation.
The small finite-volume effects are corrected using the method by Hansen and Patella [48,49], see [4] for the details of our implementation.We find the size of the correction to be at most two permille on the ensembles included in our computation.Finite-volume corrections due to kaon loops are non-negligible close to the SU(3) symmetric point, where pion and kaon masses are of similar order.We use the same procedure as for the effects from pions to estimate the corresponding correction.Both contributions and their sum are given in table V.A flat 25% systematic uncertainty is assigned to the corrections.
As explained above, we scan a variety of ansätze for the chiral-continuum extrapolation to determine (a 3,3 µ ) SD sub at the physical point.The best fit according to the AIC for the data of set 1 with the LL discretization is shown on the left-hand side of Figure 2. The minimized χ 2 is 7.64 for 11 degrees of freedom (dof).We show the dependence of (a 3,3 µ ) SD sub (25 GeV 2 ) on the light quark mass proxy Φ 2 , where the data has been corrected for deviations from Φ phys 4 and cuts to m π < 400 MeV and a < 0.09 fm have been applied.The data and the fit function are shown for the five values of the inverse gauge coupling that are included in this fit.As visible from the figure, the quark mass effects are mild but depend on the lattice spacing.The fit also includes cutoff effects proportional to a 2 Φ 4 which are not visible in the plot.
On the right-hand side of Figure 2, we show the approach to the continuum limit at physical quark masses, as determined from the scan over all fit models for both discretizations of the current and both improvement schemes.Each line corresponds to one fit and the opacity is related to the model weight of the fit according to the AIC.The relative size of the cutoff effects, about 10% in the most extreme case, is not larger than in our computation of the intermediate-distance window observable [4], thanks to some of the improvements described in section II C. For the same reason, the difference between the two discretizations of the current is small.
The two sets of improvement schemes lead to cutoff effects that have opposite signs.A similar behaviour has been observed for the case of the intermediate-distance window, where the extrapolations of both data sets agreed in the continuum limit, taking into account the curvature in the extrapolation of the set 2 data.In the case of (a 3,3 µ ) SD sub we observe that most of the fits from set 2 with a large model weight favour slightly larger values than the fits to the data of set 1.However, for some fits we observe a strong curvature at very small values of the lattice spacing such that they are compatible with the more benign extrapolations for set 1.In the current situation, where we do not have more information on the behaviour towards the continuum limit, e.g., by adding data at even finer lattice spacings, we decide to take the discrepancy between the two data sets into account as systematic uncertainty of our continuum extrapolation.
Whereas the systematic uncertainty of 0.21 amounts to about 0.6% of the value of (a 3,3 µ ) SD sub , we note that it is sub-permille with respect to a 3,3  µ .The description of the short-distance cutoff effects therefore does not seem to be a roadblock towards the target precision for the full HVP.
The qualitative behaviour of the continuum extrapolation does not vary with the choice of Q 2 .A smaller value of Q 2 , which removes more of the potentially difficult ultra short-distance region, does not lead to a better agreement of the results based on the two improvement schemes.FIG.3: Left: The ratio (b (3,3) GeV) 2 for the isovector channel, computed non-perturbatively in lattice QCD vs. using the perturbative Adler function (see eqs. ( 9), ( 16) and ( 19)).For visibility, the data is slightly displaced horizontally.Right: Stability of the result for (a 3,3 µ ) SD as a function of the virtuality Q 2 chosen in the perturbatively treated subtraction term.Choosing Q = ∞ corresponds to the entire (a 3,3 µ ) SD being computed in lattice QCD without any use of perturbation theory.

C. Combination of perturbative and non-perturbative evaluations in the isovector channel
The non-perturbative evaluation of (a hvp µ ) SD sub at several values of Q 2 gives us the possibility to monitor the convergence of perturbation theory and to directly compare the perturbative and non-perturbative evaluations of the difference in eq.(19).We choose the reference energy Q ref = 8 GeV where the convergence of perturbation theory is expected to be excellent.
In the non-perturbative evaluation, we compute the left hand side of eq. ( 19).We get compatible results when performing the difference either in the continuum limit or when subtracting ensemble by ensemble and extrapolating the difference.We choose the second version for the comparison with perturbation theory in Figure 3.It can be seen that both approaches give essentially the same results down to Q = 4 GeV, confirming the good convergence of the perturbative series.
On the right-hand side of Figure 3 we depict the final quantity (a 3,3 µ ) SD depending on the choice of Q.Any dependence on Q should be removed after combining lattice and perturbation theory and indeed all results are entirely compatible with each other.For our final estimate, we use Q = 5 GeV and obtain (a 3,3  µ ) SD = 43.06(4) stat (21) syst (3) scale [22] .
We also show a result for the direct evaluation of (a 3,3 µ ) SD , without any subtraction in the kernel function, in Figure 3, where we have ignored the possibility of cutoff effects of O(a 2 log(a)) [19] being present in the lattice data and employed the same set of continuum extrapolations as for the subtracted quantity.The agreement with our final result is an indication that the size of these effects is small for this specific observable.
FIG. 4: Illustration of fits to ∆ ls (a µ ) SD .Left: Best fit, according to its model weight, to the data for set 2, LL with χ 2 /dof = 13.07/13.The data is corrected for deviation from physical Φ 4 .The coloured lines show the evaluation of the fit function at finite lattice spacing.The black line, together with the gray band, show the dependence on Φ 2 in the continuum limit.Right: Projection onto the lattice spacing dependence and physical quark masses for the best fits in the four categories entering the model average.

D. Non-perturbative evaluation of ∆ ls (a µ ) SD
To compute the isoscalar contribution to (a hvp µ ) SD , it suffices to evaluate ∆ ls (a µ ) SD as given by eq. ( 21).The quark-disconnected loops for light and strange quarks, computed as outlined in [15], enter here and may be precisely determined in the short-distance region.By definition, the splitting ∆ ls (a µ ) SD vanishes at the SU(3) symmetric point and is expected to depend at leading order on m s − m l .We take this knowledge into account for the chiral and continuum extrapolation.
and again build a variety of ansätze by setting coefficients to zero.The ensembles with SU(3) symmetry are excluded from these fits.All cutoff effects are suppressed by Φ δ close to the SU(3) symmetric point.
On the left-hand side of figure 4 we show the Φ 2 dependence for the data and a typical fit for set 2 and the LL discretization.The constraint ∆ ls (a µ ) SD = 0 for Φ 4 = 3  2 Φ 2 at the SU(3) symmetric point is visible by the intersection of the curves which denote the Φ 2 dependence at finite lattice spacing and in the continuum, respectively.Towards the physical point, cutoff effects grow since the suppression is lifted.On the right-hand side of the figure, we show the approach to the continuum for all fits that are contained in the scan over the different models, evaluated at physical light and strange quark masses.The data from set 1 approach the continuum value from below whereas the cutoff effects for data set 2 have the opposite sign.
With the result 1 3 ∆ ls (a µ ) SD = −0.495(7) stat (34) syst (4) scale [36] , from an average over the fit models, eq. ( 21) and the result for the isovector contribution from eq. ( 43), we arrive at for the isoscalar contribution.We note that this result has significant correlation with the result for (a 3,3 µ ) SD , which is taken into account when combining both in (a hvp µ ) SD .For comparison with other lattice results, we combine the results for the isovector and isoscalar contributions in eqs.( 43) and ( 46) to compute the strange-connected contribution, 1 9 (a s,s µ ) SD = 9.072 (10) stat (58) syst (3) scale [60] .
The disconnected contribution from light and strange quark loops is found to be irrelevant with respect to the precision quoted here.A direct evaluation of the strange-connected contribution via the same strategy as for the charm-connected contribution gives a result that is fully compatible with eq. ( 47).

E. Evaluation of the charm-connected contribution
To compute the charm-connected contribution on the lattice, we start by evaluating (a c,c µ ) SD sub (Q 2 ).The tuning of the charm quark hopping parameter to match the mass of the D s meson of 1968.47MeV and the determination of the renormalization constant for the local charm current in a massive renormalization scheme have been described in [4,20].After the continuum extrapolation, we perform a small correction to adapt the tuning to the updated value of the scale setting quantity t phys 0 from [21].In contrast to the isovector contribution, we find very good agreement of the continuum extrapolations of the four data sets, despite having significantly different cutoff effects.We illustrate the cutoff effects based on the scan over the fit models for each of the four data sets on the left-hand side of Figure 5.For our final value, in line with our previous work, we use only the LC discretization of the current, since it exhibits significantly smaller cutoff effects.
The effect of massive quarks in the contribution b (c,c) may be computed on the lattice without potentially dangerous log-enhanced cutoff effects.The size of this contribution strongly depends on the choice of Q and at large values of Q, the observable is dominated by distances < 0.3 fm.On the right-hand side of Figure 5 we show the approaches to the continuum, based on the scan of the fit models.We evaluate ∆ lc b(Q 2 ) for the same values of ) computed non-perturbatively in lattice QCD vs. using massive perturbation theory.For visibility, the data is slightly displaced horizontally.Right: Stability of the result for (a c,c µ ) SD as a function of the virtuality Q 2 chosen in the subtraction term.Choosing Q = ∞ corresponds to the entire (a c,c µ ) SD being computed in lattice QCD without any use of perturbation theory.
Q 2 as (a 3,3 µ ) SD to explore the parameter space and to be able to compare to the evaluation of the same observable using massive perturbation theory, see table II.
The perturbative evaluation is based on the expansion in m 2 c /Q 2 of the vacuum polarization, computed to O(α 2 s ) in [50].We have kept expansion terms up to O((m 2 c /Q 2 ) 4 ) in the perturbative orders α 0 s and α s , while we kept expansion terms up to O((m , we evaluate α s as well as the MS charm mass at a scale µ 2 given by the geometric mean of Q 2 and Q 2 /4, using the FLAG'21 [46] value of Λ MS and starting from the PDG value m c (µ = m c ) = 1.27GeV [51].The number indicated in brackets in the b (c,c) column of table II represents the contribution of the O(α 2 s ) term as a way to indicate the uncertainty of the prediction.The quark-mass dependent contribution to b (c,c) (Q 2 ) is compared to lattice data in Fig. 6.
On the left-hand side of Figure 6, we show the comparison of the perturbative prediction and the lattice result for ∆ lc b(Q 2 ) in the form proposed in eq.(19).We choose Q ref = ∞ where ∆ lc has to vanish.It is visible that at small Q, despite giving similar central values, the uncertainty on the perturbative result is significantly larger than the uncertainty on the lattice evaluation, which is dominated by the systematic uncertainty from the continuum extrapolation.We note that the relative size of this uncertainty grows with increasing Q, whereas its absolute size decreases.
For our final result, we combine the non-perturbative evaluations of (a c,c µ ) SD sub (Q 2 ) and ∆ lc b(Q 2 ) with the perturbative result for b (3,3) (Q 2 ).The combination for several values of Q 2 is shown on the right-hand side of Figure 6.The residual dependence on Q is negligible with respect to the uncertainties.For our final result, we again choose Q = 5 GeV and with ∆ lc b(25 GeV 2 ) = −2.42(4)stat (10) syst (3) scale [11] , we arrive at 4 9 (a c,c µ ) SD = 11.53 (13) stat (23) syst (11) scale [30] . .The coloured lines show evaluations of the fit function at finite lattice spacing and the grey area show the result in the continuum limit.Left: Fit to 4  9 (a c,c µ ) SD disc with χ 2 /dof = 17.60/16.Right: Fit to 2 3 √ 3 (a c,8 µ ) SD disc with χ 2 /dof = 9.20/12.

F. The charm-disconnected contributions
The quark-disconnected contributions including charm quarks have not been included in our previous results for a hvp µ [20] and (a hvp µ ) ID [4] because they are expected to be numerically very small and mostly contained in the short-distance region.In this work, we compute the contributions G (c,8) and G (c,c) disc in a partially quenched setup, see Appendix C of [15] for an explanation of the computational setup.The charm quark hopping parameter is set to the same value as in the quark-connected case.To avoid mixing with the singlet current at O(a) in the case of G (c,8) , and anticipating smaller cutoff effects in the case of a heavy quark (see section III E), we evaluate the correlation functions with conserved vector currents at source and sink.
We find both contributions to be very small and, at finite lattice spacing, non-zero only below a distance of about 1 fm.Given the smallness of these contributions, we perform the analysis without subtractions in the kernel function, ignoring the possibility of log-enhanced cutoff effects.
The charm-charm contribution is dominated by cutoff effects, already visible by comparing different discretizations of the current on a single ensemble.We find significantly different integrands and even a change in sign when employing the LL formulation of the current instead of CC or LC.Accordingly, the relative size of the cutoff effects is found to be large in the continuum extrapolation.On the left-hand side of Figure 7 we show an exemplary chiral-continuum extrapolation for this contribution.As can bee seen in the figure, the sign changes from positive to negative when taking the continuum limit.We attribute the slight light quark mass dependence to the matching of the charm quark via the mass of the D s meson and the variation of the strange quark mass in our set of ensembles.At the physical point we find after averaging over the fit models the value 4 9 (a c,c µ ) SD disc = −0.0010(18) stat (32) syst (1) scale [37] , and thus only an upper limit on the magnitude of this contribution.
Since the contribution (a c,8 µ ) SD is SU(3) suppressed, we employ fits according to eq. ( 44).An example fit is shown on the right-hand side of Figure 7. Cutoff effects are suppressed close to the SU(3) symmetric point but significant at physical quark mass.The chiral-continuum extrapolated value is compatible with zero.We obtain µ ) SD disc = 0.0020 (17) stat (25) syst (0) scale [31] , from averaging over all fit models.Both contributions are thus negligible compared to (a hvp µ ) SD (and our final uncertainty).Naturally, this carries over to a hvp µ .Summing the two results above and adding the disconnected contribution that enters our final result via (a 8,8  µ ) SD , we quote (a hvp µ ) SD disc = 0.0013 (25) stat (41) syst (5) scale [49] to allow for a comparison with other works.

G. Isospin-breaking corrections
We apply two complementary methods to determine isospin-breaking effects in (a hvp µ ) SD .The first employs massless QCD perturbation theory to compute the leading QED effects.The second approach is based on the explicit calculation of isospin-breaking corrections due to unequal up and down quark masses and electric charges, following the same method as in our previous publication on the intermediate-distance window observable [4].
To compute the leading QED effects at short distances in massless QCD perturbation theory, we start from the expression for the relative correction to the Adler function (see [52], where the correction is given for the R ratio) Note that in this notation D QCD (Q 2 ) = 3(1+(α s /π)+O(α 2 s )) f =u,d,s Q 2 f is the pure massless QCD expression, while D QCD,1 internal γ (Q 2 ) is the contribution from massless QCD containing exactly one internal photon line.The ratio (54) amounts to 0.51 × 10 −3 for α s = 0.30 for the (u, d, s) contribution, and to 0.57 × 10 −3 if one restricts one's attention to the (u, d) sector.These are very small corrections indeed, yielding the estimate [(a 3,3  µ ) SD + 1 3 (a 8,8 µ ) SD ] 1 internal γ = 0.03.For the QED correction to the charm contribution, we use the KS spectral function expressed in terms of the MS charm mass.This means that the latter mass is kept fixed as the electromagnetic interaction is 'turned on'.The free-quark level contribution of the charm is 4  9 (a c,c µ ) SD = 11.6 using the mass of m c (m c ) = 1.27GeV, while the relative QED correction to that is 0.83 × 10 −3 .In absolute terms, this means For our second approach we have computed (a hvp µ ) SD in QCD+QED on a subset of our ensembles using Monte Carlo reweighting combined with the leading-order perturbative expansion of QCD+QED about the isosymmetric theory [53][54][55][56][57].The same method was applied in Ref. [4].In order to determine the dependence of isospin-breaking corrections to the short-distance window observable on the pion mass and the lattice spacing, we have used eight ensembles (A654, H102, N101, N452, N451, D450, N203 and N200), covering a much wider range of pion masses (220−350 MeV) and lattice spacings (0.064−0.097 fm) compared to our previous calculation [4].As before, we define our hadronic renormalization scheme in terms of m 2 π 0 , m 2 π 0 and the fine-structure constant α.Neglecting isospin-breaking effects in the lattice scale and in the quark sea, as well as quark-disconnected contributions in the calculation of the relevant correlation functions, we realize our scheme by matching the values of m 2 π 0 and m 2 K + + m 2 K 0 − m 2 π + in QCD+QED to those in the isosymmetric theory and by setting m 2 π 0 to its experimental value.Further technical details are described in section VI of Ref. [4].
Our results for the relative size of the isospin-breaking corrections to the connected light and strange contributions to (a hvp µ ) SD are shown in the upper panel of Figure 8 for the conserved-local (left) and local-local (right) combinations of vector currents.We make two observations: first, the relative corrections are very small and typically amount to a few per-mil.Second, there is no clear dependence of the results on the pion mass and the lattice spacing, which would allow for a systematic extrapolation to the physical point.On the other hand, our results for the meson mass splittings m K + − m K 0 and m π + − m π 0 , shown in the lower panels of Figure 8, and which include an additional ensemble at almost physical pion mass (D452), show a consistent trend towards the corresponding experimental values.We are, therefore, confident that the typical size of isospin-breaking corrections to (a hvp µ ) SD has been determined correctly.Taking into account the spread in the results and allowing for a generous error, we represent the data shown in the upper panels of Figure 8 by (1.5 ± 1.5) × 10 −3 .We can then work out the absolute correction to the connected light and strange quark contribution, by multiplying the sum of eqs.( 43) and (47) with 0.0015, which yields [(a 3,3 µ ) SD + 1 9 (a s,s µ ) SD ] × 0.0015 = 0.085 .
Adding the QED correction to the charm quark contribution of 0.01 from eq. ( 55) and assuming an uncertainty of 100% for the latter, we arrive at which we quote as the total isospin-breaking correction to (a hvp µ ) SD .This number is larger but of the same order of magnitude compared to the QED correction to [(a 3,3  µ ) SD + 1 3 (a 8,8 µ ) SD + 4 9 (a c,c µ ) SD ] estimated in massless QCD perturbation theory at leading order and denoted by the light blue crosses in Fig. 8, which further corroborates our findings.

H. Further, small contributions
In this section, we present our (mostly perturbative) estimates for two small contributions that we did not evaluate directly in lattice QCD.
For the contribution of the b quark, (a b,b µ ) SD , the perturbative spectral function is known exactly to O(α s ) from the calculation of Källén-Sabry (KS) [58], described for instance in [59].The O(α 2 s ) corrections were calculated in [60].We have evaluated the KS prediction, improved by substituting the pole quark-mass appearing in R(s) by its MS counterpart and evaluating the latter at the scale s.Using the PDG bottom mass of 4.18 GeV, this results in the estimate 1  9 (a b,b µ ) SD = 0.31.We also note the NRQCD based lattice calculation [61], which finds 1  9 a b,b µ = 0.271(37), while the phenomenological estimate [62] based on sum rules obtains 0.30 (2).We have checked that the bottom contribution to 1  9 (a b,b µ − (a b,b µ ) SD ) is on the order of 0.003.We thus arrive at our estimate of As for the effect of the missing charm sea quarks in our lattice calculation, we estimate its order of magnitude based on perturbation theory, rather than on D meson loops as we previously did for the intermediate window in [4], since (a hvp µ ) SD involves relatively short distances.Using the perturbative results of [63], we estimate the effect to be at the level of [(a 3,3  µ ) SD + 1 3 (a 8,8 µ ) SD ] pert.charm sea quark effect = 0.02.
Note that the perturbative estimate of the connected charm contribution is in rather good agreement with the lattice calculation (see the text above eq.( 55)).Nevertheless, since eq.( 59) represents only the leading perturbative prediction for the effect under scrutiny, we will conservatively assign an uncertainty of ∆[(a 3,3 µ ) SD + 1 3 (a 8,8 µ ) SD ] charm quenching = 0.10 to our calculation due to the quenching of the charm quark.

IV. CONCLUSION
Combining the results of sections III, we arrive at the result (a hvp µ ) SD = 68.85(14) stat (39) syst (15) scale [45] , where all contributions and sources of uncertainty have been taken into account.As outlined in the main text, we have evaluated the subtracted quantity of eq. ( 17) for Q 2 = 25 GeV 2 .For this choice, 26% of the result in eq. ( 61) are due to massless perturbation theory applied in the spacelike domain.On the left-hand side of Figure 9 we show results for (a hvp µ ) SD for several choices of the virtuality, the grey band denoting our final estimate.No residual dependence on the choice for Q can be resolved.At Q = 4 GeV 41% of the result would stem from massless perturbation theory, as given by Table II, whereas the contribution at Q = 8 GeV would only be 10%.We find a noticeable, but not significant difference of our results using (a hvp µ ) SD sub with respect to a direct evaluation of (a hvp µ ) SD , as indicated by the leftmost data point in the figure .We have shown that, in our data set, (a hvp µ ) SD depends much more strongly on the lattice spacing than on the quark masses (see Fig. 2).The very fine lattices employed here were thus instrumental in controlling the continuum limit, even though most of them correspond to heavier-than-physical pion masses.The final uncertainty is dominated by systematic uncertainties, mostly from the variation of models selected to perform the continuum extrapolation.On the right-hand side of Fig. 9, we show the contributions to the squared uncertainty.The contribution labelled 'statistics' corresponds to the statistical uncertainty of 0.14 in eq. ( 61).The dependence on the scale setting quantity, and therefore also the resulting contribution to the final uncertainty, is dominated by the charm quark contribution and contains the matching with the mass of the D s meson.
In table III we collect the dependence of specific (intermediate) results with respect to the quantities that define our scheme of isospin symmetric QCD.The information can be used to adapt the results to a different scheme, provided that the differences in the input quantities are small.For example, a shift in the scale setting quantity t phys 0 to the current N f = 2 + 1 + 1 FLAG average with the central value (t 0 ) phys N f =2+1+1 = 0.14292 fm would result in the adapted value (a hvp µ ) SD = 69.00(45).FIG.10: Comparison of our results for the short-distance window observable (a hvp µ ) SD (right panel) and the light-quark connected contribution (left panel) to other recent lattice calculations by ETM [5] and RBC/UKQCD [6].Our results are denoted by green circles and band, while the data-driven estimate of Colangelo et al. [7] is shown as the red circle and band.
A compilation of recent results for the short-distance window observable [5][6][7], including the estimates from this work, is shown in Figure 10.Our result (61) is in good agreement with the dispersive estimate from [7], (a hvp µ ) SD = 68.4(5).Our result is also in good agreement with the lattice result from the ETM collaboration [5], (a hvp µ ) SD = 69.27(34).In the introduction, a schematic scenario was discussed in which the R ratio corresponding to lattice calculations is enhanced by six percent relative to the e + e − data entering the determination of [7] in the interval 600 MeV < √ s < 900 MeV.Such an enhancement would increase (a hvp µ ) SD by about one unit.Both the ETM and our result are also consistent with this scenario, which thus remains a working hypothesis worth challenging.
Our result for the light-quark connected contribution, 10 9 (a 3,3 µ ) SD = 47.84 (24) is consistent with the ETM result [5] 48.24 (20) and the RBC/UKQCD result [6] 48.70(52)(59) at their favoured isospin-symmetric point.The dependence of (a hvp µ ) SD on the precise choice of this point is currently smaller than the statistical errors.
Finally, the subtraction strategy for the short-distance effects adopted here can be reused for the full a hvp µ (the analogue of eq. ( 15) being a hvp µ = (a hvp µ ) sub (Q 2 ) + b (γ,γ) (Q 2 )), as well as for other observables related to the vacuum polarization, such as the running of the electromagnetic coupling.The appeal of this subtraction is that it removes the leading, O(x 4 0 ) term from the kernel, while at the same being safe to compute in perturbation theory, since it only involves spacelike photon virtualities on the order of Q 2 .TABLE V: Finite-size effects for (a hvp µ ) SD in the isovector channel, in units of 10 −10 .The correction is obtained using the Hansen-Patella method and split according to eq. (15).The first part enters our calculation in the non-perturbative determination of (a 3,3 µ ) SD sub (25 GeV 2 ), the second one in the determination of ∆ lc b(25 GeV 2 ).The columns contain the correction due to pions and kaons wrapping around the torus.On ensembles with SU(3) symmetry, the correction is contained in the correction due to pions.The column labelled 'total' gives the sum of the two contributions.When including the finite-volume corrections in our analysis, we assign a flat 25% relative uncertainty to the correction.

FIG. 2 :
FIG.2: Illustration of fits to (a3,3  µ ) SD sub (25 GeV 2 ).Left: Best fit, according to its model weight, to the data for set 1, LL.The data is corrected for deviation from physical Φ 4 .The coloured lines show the evaluation of the fit function at finite lattice spacing.The black line, together with the grey band, show the dependence on Φ 2 in the continuum limit.Right: Approaches to the continuum limit for four sets of data based on the improvement schemes of set 1 and 2 and the LL and LC discretizations of the current based on a scan over fit models for (a 3,3 µ ) SD sub (25 GeV 2 ).Each line shows the result from one single fit and the opacity of the lines corresponds to the weight of the fit in the model average.Dashed vertical lines indicate the lattice spacings used in this work.

FIG. 6 :
FIG.6: Left: The difference −4  9 Q 2 ∆ lc b(Q 2 ) computed non-perturbatively in lattice QCD vs. using massive perturbation theory.For visibility, the data is slightly displaced horizontally.Right: Stability of the result for (a c,c µ ) SD as a function of the virtuality Q 2 chosen in the subtraction term.Choosing Q = ∞ corresponds to the entire (a c,c µ ) SD being computed in lattice QCD without any use of perturbation theory.

FIG. 7 :
FIG.7: Illustration of typical fits to the charm-disconnected contributions.The data is corrected for deviations from Φ phys 4

55 FIG. 8 :
FIG.8:The relative strong and electromagnetic isospin-breaking corrections to the connected light and strange quark contribution computed using the conserved-local (upper left panel) and local-local (upper right panel) discretizations of the vector current.The light blue crosses denote the prediction from massless QCD (see the main text).The dashed line and the grey area show our final estimate.The two panels at the bottom show the mass splittings of charged and neutral kaons and pions, respectively.The black stars represent the experimental values.

FIG. 9 :
FIG.9: Left: Residual Q dependence of (a hvp µ ) SD .The grey area shows the final result.Right: Contributions to the variance of (a hvp µ ) SD as evaluated for Q = 5 GeV.

TABLE I :
Parameters of the simulations: the bare coupling β = 6/g 2 0 , the temporal boundary conditions, open (o) or anti-periodic (p), the lattice dimensions, the lattice spacing a in physical units based on

TABLE II :
Lattice spacing dependence from a scan over fit models for four data sets.The lattice spacings used in this work are indicated by vertical lines.Left: fits to (a c,c µ ) SD sub (36 GeV 2 ).Right: fits to − 4 9 ∆ lc b(25 GeV 2 ).

TABLE III :
Dimensionless scheme dependencies of observable O with respect to the quantity S according to S ∂O ∂S .The four quantities S define the scheme of isoQCD in this work.Their central values are √ t 0 = 0.1443 fm, m π = 134.9768MeV, m K = 495.011MeV, m Ds = 1968.47MeV.

TABLE VII :
Values of the ∆ ls (a µ ) SD contribution in units of 10 −10 , for the local-local (LL) and for the local-conserved (CL) discretizations of the correlation function, as described in the main text.The finite-size correction has been applied.− 1 3 ∆ ls (a µ ) SD -Set 1 − 1 3 ∆ ls (a µ ) SD -Set 2 id

TABLE VIII :
Values of the ∆ lc b contribution with Q 2 = 25 GeV 2 in units of 10 −10 , for the locallocal (LL) and for the local-conserved (CL) discretizations of the correlation function, as described in the main text.The finite-size correction has been applied.