Window observable for the hadronic vacuum polarization contribution to the muon $g-2$ from lattice QCD

Euclidean time windows in the integral representation of the hadronic vacuum polarization contribution to the muon $g-2$ serve to test the consistency of lattice calculations and may help in tracing the origins of a potential tension between lattice and data-driven evaluations. In this paper, we present results for the intermediate time window observable computed using O($a$) improved Wilson fermions at six values of the lattice spacings below 0.1\,fm and pion masses down to the physical value. Using two different sets of improvement coefficients in the definitions of the local and conserved vector currents, we perform a detailed scaling study which results in a fully controlled extrapolation to the continuum limit without any additional treatment of the data, except for the inclusion of finite-volume corrections. To determine the latter, we use a combination of the method of Hansen and Patella and the Meyer-Lellouch-L\"uscher procedure employing the Gounaris-Sakurai parameterization for the pion form factor. We correct our results for isospin-breaking effects via the perturbative expansion of QCD+QED around the isosymmetric theory. Our result at the physical point is $a_\mu^{\mathrm{win}}=(237.30\pm0.79_{\rm stat}\pm1.22_{\rm syst})\times10^{-10}$, where the systematic error includes an estimate of the uncertainty due to the quenched charm quark in our calculation. Our result displays a tension of 3.9$\sigma$ with a recent evaluation of $a_\mu^{\mathrm{win}}$ based on the data-driven method.


I. INTRODUCTION
The anomalous magnetic moment of the muon, a µ , plays a central role in precision tests of the Standard Model (SM).The recently published result of the direct measurement of a µ by the Muon g − 2 Collaboration [1] has confirmed the earlier determination by the E821 experiment at BNL [2].When confronted with the theoretical estimate published in the 2020 White Paper [3], the combination of the two direct measurements increases the tension with the SM to 4.2σ.The SM prediction of Ref. [3] is based on the estimate of the leading-order hadronic vacuum polarization (HVP) contribution, a hvp µ , evaluated from a dispersion integral involving hadronic cross section data ("data-driven approach") [4][5][6][7][8][9], which yields a hvp µ = (693.1 ± 4.0) × 10 −10 [3].The quoted error of 0.6% is subject to experimental uncertainties associated with measured cross section data.
Lattice QCD calculations for a hvp µ [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] as well as for the hadronic light-by-light scattering contribution a hlbl µ [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39] have become increasingly precise in recent years (see [40][41][42] for recent reviews).Although these calculations do not rely on the use of experimental data, they face numerous technical challenges that must be brought under control if one aims for a total error that can rival or even surpass that of the data-driven approach.In spite of the technical difficulties, a first calculation of a hvp µ with a precision of 0.8% has been published recently by the BMW collaboration [20].Their result of a hvp µ = (707.5 ± 5.5) × 10 −10 is in slight tension (2.1σ) with the White Paper estimate and reduces the tension with the combined measurement from E989 and E821 to just 1.5σ.This has triggered several investigations that study the question whether the SM can accommodate a higher value for a hvp µ without being in conflict with low-energy hadronic cross section data [43] or other constraints, such as global electroweak fits [44][45][46][47].At the same time, the consistency among lattice QCD calculations is being scrutinized with a focus on whether systematic effects such as discretization errors or finite-volume effects are sufficiently well controlled.Moreover, when comparing lattice results for a hvp µ from different collaborations, one has to make sure that they refer to the same hadronic renormalization scheme that expresses the bare quark masses and the coupling in terms of measured hadronic observables.
Given the importance of the subject and in view of the enormous effort required to produce a result for a hvp µ at the desired level of precision, it has been proposed to perform consistency checks among different lattice calculations in terms of suitable benchmark quantities that suppress, respectively enhance individual systematic effects.These quantities are commonly referred to as "window observables", whose definition is given in section II.
In this paper we report our results for the so-called "intermediate" window observables, for which the short-distance as well as the long-distance contributions in the integral representation of a hvp µ are reduced.This allows for a straightforward and highly precise comparison with the results from other lattice calculations and the data-driven approach.This constitutes a first step towards a deeper analysis of a possible deviation between lattice and phenomenology.Indeed, our findings present further evidence for a strong tension between lattice calculations and the data-driven method.At the physical point we obtain a win µ = (237.30± 1.46) × 10 −10 (see Eq. ( 45)) for a detailed error budget), which is 3.9σ above the recent phenomenological evaluation of (229.4 ± 1.4) × 10 −10 quoted in Ref. [48].
This paper is organized as follows: We motivate and define the window observables in Sect.II, before describing the details of our lattice calculation in Sect.III.In Sect.IV we discuss extensively the extrapolation to the physical point, focussing specifically on the scaling behavior, and present our results for different isospin components and the quark-disconnected contribution.Sections V and VI describe our determinations of the charm quark contribution and of isospin-breaking corrections, respectively.Our final results are presented and compared to other determinations in Sect.VII.In-depth descriptions of technical details and procedures, as well as data tables, are relegated to several appendices.Details on how we correct for mistunings of the chiral trajectory are described in Appendices A and B, the determination of finite-volume corrections is discussed in Appendix C, while the estimation of the systematic uncertainty related to the quenching of the charm quark is presented in Appendix D. Ancillary calculations of pseudoscalar masses and decay constants that enter the analysis are described in Appendix E. Finally, Appendix F contains extensive tables of our raw data.

II. WINDOW OBSERVABLES
The most widely used approach to determine the leading HVP contribution a hvp µ in lattice QCD is the "time-momentum representation" (TMR) [49], i.e.
where G(t) is the spatially summed correlation function of the electromagnetic current K(t) is a known kernel function (see Appendix B of Ref. [10]), and the integration is performed over the Euclidean time variable t.By considering the contributions from the light (u, d), strange and charm quarks to G(t) one can perform a decomposition of a hvp µ in terms of individual quark flavors.It is also convenient to consider the decomposition of the electromagnetic current into an isovector (I = 1) and an isoscalar (I = 0) component according to where the ellipsis in the first line denotes the missing charm and bottom contributions.
One of the challenges in the evaluation of a hvp µ is associated with the long-distance regime of the vector correlator G(t).Owing to the properties of the kernel K(t), the integrand K(t)G(t) has a slowly decaying tail that makes a sizeable contribution to a hvp µ in the region t > ∼ 2 fm.However, the statistical error in the calculation of G(t) increases exponentially with t, which makes an accurate determination a difficult task.Furthermore, it is the long-distance regime of the vector correlator that is mostly affected by finite-size effects.
The opposite end of the integration interval, i.e. the interval t < ∼ 0.4 fm, is particularly sensitive to discretization effects which must be removed through a careful extrapolation to the continuum limit, possibly involving an ansatz that includes sub-leading lattice artefacts, especially if one is striving for sub-percent precision.
At this point it becomes clear that lattice results for a hvp µ are least affected by systematic effects in an intermediate subinterval of the integration in Eq. ( 1), as already recognized in [49].This led the authors of Ref. [13] to introduce three "window observables", each defined in terms of complementary sub-domains with the help of smoothed step functions.To be specific, the short-distance (SD), intermediate distance (ID) and long-distance (LD) window observables are given by where ∆ denotes the width of the smoothed step function Θ defined by Θ(t, t , ∆) The widely used choice of intervals and smoothing width that we will follow is t 0 = 0.4 fm, t 1 = 1.0 fm and ∆ = 0.15 fm.(8) The original motivation for introducing the window observables in Ref. [13] was based on the observation that the relative strengths and weaknesses of the lattice QCD and the R-ratio approach complement each other when the evaluations using either method are restricted to nonoverlapping windows, thus achieving a higher overall precision from their combination.Since then it has been realized that the window observables serve as ideal benchmark quantities for assessing the consistency of lattice calculations, since the choice of sub-interval can be regarded as a filter for different systematic effects.Furthermore, the results can be confronted with the corresponding estimate using the data-driven approach.This allows for high-precision consistency checks among different lattice calculations and between lattice QCD and phenomenology.
In this paper, we focus on the intermediate window and use the simplified notation We remark that the observable a win µ , which accounts for about one third of the total a hvp µ , can be obtained from experimental data for the ratio via the dispersive representation of the correlator (2) [49].How different intervals of centerof-mass energy contribute to the different window observables in the data-driven approach is investigated in Appendix B; similar observations have already been made in Refs.[48,50,51].
For the intermediate window a win µ , the relative contribution of the region √ s < 600 MeV is significantly suppressed as compared to the quantity a hvp µ .Instead, the relative contribution of the region √ s > 900 MeV, including the φ meson contribution, is somewhat enhanced 1 .
Interestingly, the region of the ρ and ω mesons between 600 and 900 MeV makes about the same fractional contribution to a win µ as to a hvp µ , namely 55 to 60%.Thus if the spectral function associated with the lattice correlator G(t) was for some reason enhanced by a constant factor (1 + ) in the interval 600 < √ s/MeV < 900 relative to the experimentally measured spectral function R(s)/(12π 2 ), it would approximately lead to an enhancement by a factor (1 + 0.6 ) of both a hvp µ and a win µ .Finally, we note that the relative contributions of the three √ s intervals are rather similar for a win µ as for the running of the electromagnetic coupling from III. CALCULATION OF a win µ ON THE LATTICE A. Gauge ensembles Our calculation employs a set of 24 gauge ensembles generated as part of the CLS (Coordinated Lattice Simulations) initiative using N f = 2 + 1 dynamical flavors of non-perturbatively O(a) improved Wilson quarks and the tree-level O(a 2 ) improved Lüscher-Weisz gauge action [52].The gauge ensembles used in this work were generated for constant average bare quark mass such that the improved bare coupling g0 [53] is kept constant along the chiral trajectory.Six of the ensembles listed in Table I Pion masses lie in the range m π ≈ 130 − 420 MeV.Seven of the ensembles used have periodic (anti-periodic for fermions) boundary conditions in time, while the others admit open boundary conditions [54].All ensembles included in the final analysis satisfy m π L > ∼ 4. Finite-size effects can be checked explicitly for m π = 280 and 420 MeV, where in each case two ensembles with different volumes but otherwise identical parameters are available.The ensembles with volumes deemed to be too small are marked by an asterisk in Table I and are excluded from the final analysis.
The QCD expectation values are obtained from the CLS ensembles by including appropriate reweighting factors, including a potential sign of the latter [55].A negative reweighting factor, which originates from the handling of the strange quark, is found on fewer than 0.5% of the gauge field configurations employed in this work. 1Contributions as massive as the J/ψ, however, make again a smaller relative contribution to a win µ than to a hvp µ .
TABLE I: Parameters of the simulations: the bare coupling β = 6/g 2 0 , the lattice dimensions, the lattice spacing a in physical units extracted from [56], the pion and kaon masses and the physical size of the lattice, the number of gauge field configurations used for the connected light-and strange-quark contributions (penultimate column) and for the disconnected contribution (last column).Ensembles with an asterisk are not included in the final analysis but used to control finite-size effects.The ensembles A653, A654, B450, N451, D450, D452, and E250 have periodic boundary conditions in time, all others have open boundary conditions.

Id
β For the bulk of our pion masses, down to the physical value, results were obtained at four values of the lattice spacing in the range a = 0.050 − 0.086 fm.At and close to the SU(3) fsymmetric point, four more ensembles have been added that significantly extend the range of available lattice spacings to a = 0.039 − 0.099 fm, which allows us to perform a scaling test with unprecedented precision.

B. Renormalization and O(a)-improvement
To reduce discretization effects, on-shell O(a)-improvement has been fully implemented.CLS simulations are performed using a non-perturbatively O(a) improved Wilson action [57], therefore we focus here on the improvement of the vector current in the (u, d, s) quark sector.To further constrain the continuum extrapolation and explicitly check our ability to remove leading lattice artefacts, two discretizations of the vector current are used, the local (L) and the point-split (C) currents where ψ denotes a vector in flavor space, λ are the Gell-Mann matrices, and U µ (x) is the gauge link in the direction μ associated with site x.With the local tensor current defined as , the improved vector currents are given by where ∂ is the symmetric discrete derivative ∂ν f V have been determined non-perturbatively in Ref. [58] by imposing Ward identities in large volume ensembles and independently in [59] using the Schrödinger functional (SF) setup.The availability of two independent sets allows us to perform detailed scaling tests, which is a crucial ingredient for a fully controlled continuum extrapolation.
The conserved vector current does not need to be further renormalized.For the local vector current, the renormalization pattern, including O(a)-improvement, has been derived in Ref. [60].Following the notations of Ref. [58], the renormalized isovector and isoscalar parts of the electromagnetic current read where J 0 µ = 1 2 ψγ µ ψ is the flavor-singlet current and Here, m q,l and m q,s are the subtracted bare quark masses of the light and strange-quarks respectively defined in Appendix E and m av q = (2m q,l + m q,s )/3 stands for the average bare quark mass.The renormalization constant in the chiral limit, Z V , and the improvement coefficients b V and b V , have been determined non-perturbatively in Ref. [58].Again, independent determinations using the SF setup are available in [59,61].The coefficient f V , which starts at order g 6 0 in perturbation theory [58], is unknown but expected to be very small and is therefore neglected in our analysis.
Thus, in addition to having two discretizations of the vector current, we also have at our disposal two sets of improvement coefficients that can be used to benchmark our continuum extrapolation: • Set 1 : using the improvement coefficients obtained in large-volume simulations in Ref. [58].
• Set 2 : using Z V and c V from [59], b V and b V from [61], using the SF setup.Note, in particular, that the improvement coefficients c V , b V and b V have an intrinsic ambiguity of order O(a).Thus, for a physical observable, we expect different lattice artefacts at order O(a n ) with n ≥ 2. This will be considered in Section IV C.

C. Correlation functions
The vector two-point correlation function is computed with the local vector current at the source and either the local or the point-split vector current at the sink.The corresponding renormalized correlators are with the improved correlators In the absence of QED and strong isospin breaking, there are only two sets of Wick contractions, corresponding to the quark-connected part and the quark-disconnected part of the vector two-point functions.The method used to compute the connected contribution has been presented previously in [17].In this work we have added several new ensembles and have significantly increased our statistics, especially for our most chiral ensembles.The method used to compute the disconnected contribution involving light and strange quarks is presented in detail in Ref. [62].Note that we neglect the charm quark contribution to disconnected diagrams in the present calculation.

D. Treatment of statistical errors and autocorrelations
Statistical errors are estimated using the Jackknife procedure with blocking to reduce the size of auto-correlations.In practice, the same number of 100 Jackknife samples is used for all ensembles to simplify the error propagation.In a fit, samples from different ensembles are then easily matched.
Our analysis makes use of the pion and kaon masses, their decay constants, the Wilson flow observable t 0 , as well as the Gounaris-Sakurai parameters entering the estimate of finitesize effects.These observables are always estimated on identical sets of gauge configurations and using the same blocking procedure, such that correlations are easily propagated using the Jackknife procedure.
The light and strange-quark contributions have been computed on the same set of gauge configurations, except for A654 where only the connected strange-quark contribution has been calculated.The quark-disconnected contribution is also obtained on the same set of configurations for most ensembles (see Table I).When it is not, correlations are not fully propagated; this is expected to have a very small impact on the error, since the disconnected contribution has a much larger relative statistical error.
The charm quark contribution, which is at the one-percent level, is obtained using a smaller subset of gauge configurations.Since its dependence on the ratio of pion mass to decay constant (m π /f π ) is rather flat, the error of this ratio is neglected in the chiral extrapolation of the charm contribution.
In order to test the validity of our treatment of statistical errors, we have performed an independent check of the entire analysis using the Γ-method [63] for the estimation of autocorrelation times and statistical uncertainties.The propagation of errors is based on a first-order Taylor series expansion with derivatives obtained from automatic differentiation [64].Correlations of observables based on overlapping subsets of configurations are fully propagated and the results confirm the assumptions made above.For the intermediate window observable, the contribution from the noisy tail of the correlation function is exponentially suppressed and the lattice data are statistically very precise.Thus, on each ensemble, a win µ is obtained using Eq. ( 5) after replacing the integral by a discrete sum over timeslices.Since the time extent of our correlator is far longer than t 1 = 1.0 fm, we can safely replace the upper bound of Eq. ( 5) by T /2, with T the time extent of the lattice.The results for individual ensembles are summarized in Tables VIII, IX and X.On ensemble E250, corresponding to a pion mass of 130 MeV, we reach a relative statistical precision of about two permille for both the isovector and isoscalar contributions.The integrands used to obtain a win µ are displayed in Fig. 1.
Our simulations are performed in boxes of finite volume L 3 with m π L > ∼ 4, and corrections due to finite-size effects (FSE) are added to each ensemble individually prior to any continuum and chiral extrapolation.This is the only correction applied to the raw lattice data.FSE are dominated by the ππ channel and mostly affect the isovector correlator at large Euclidean times.For the intermediate window observable, they are highly suppressed compared to the full hadronic vacuum polarization contribution.Despite this suppression, FSE in the isovector channel are not negligible and require a careful treatment.They are of the same order of FIG. 1: Integrands used to compute the intermediate window a win µ for the isovector, isoscalar and charm quark contributions.The isoscalar contribution does not include the charm quark contribution.The data has been obtained on ensemble E250, which has close-to-physical quark masses, using two local vector currents and set 1 of renormalization and improvement coefficients.magnitude as the statistical precision for our most chiral ensemble and enhanced at larger pion masses.In the isoscalar channel, FSE are included only at the SU(3) f point where m π = m K .The methodology is presented in Appendix C, and the corrections we have applied to the lattice data are given in the last column of Tables VI and V respectively for Strategy 1 and 2. In our analysis, we have conservatively assigned an uncertainty of 25% to these finite-size corrections, in order to account for any potential effect not covered by the theoretical approaches described in Appendix C. In addition to the ensembles H105 and H200 that are only used to cross-check the FSE estimate, ensembles S400 and N302 are also affected by large finite-volume corrections.We exclude those ensembles in the isovector channel.

IV. EXTRAPOLATION TO THE PHYSICAL POINT
A. Definition of the physical point in iso-symmetric QCD Our gauge ensembles have been generated in the isospin limit of QCD with m l ≡ m u = m d , neglecting strong isospin-breaking effects and QED corrections.Naively, those effects are expected to be of order O((m d − m u )/Λ QCD ) ≈ 1% and O(α) ≈ 1%, and are not entirely negligible at our level of precision.In Ref. [65], although the authors used a different scheme to define their iso-symmetric setup, those corrections have been found to be of the order of 0.4% for this window observable.A similar conclusion was reached in Ref. [13] although only a subset of the diagrams was considered.This correction will be discussed in Section VI.Only in full QCD+QED is the precise value of the observable unambiguously defined: the separation between its iso-symmetric value and the isospin-breaking correction is scheme dependent.In Section IV D, we provide the necessary information to translate our result into a different scheme.
Throughout our calculation, we define the 'physical' point in the (m π , m K ) plane by imposing the conditions [66][67][68] Inserting the PDG values [69] on the right-hand side, our physical iso-symmetric theory is thus defined by the values m π = 134.9768(5) MeV , m K = 495.011(10) MeV .
We note that since our gauge ensembles have been generated at constant sum of the bare quark masses, the linear combination (m 2 K +m 2 π /2) is approximately constant.Two different strategies are used to extrapolate the lattice data to the physical point.
a. Strategy 1 We use the gradient flow observable t 0 [70] as an intermediate scale and the dimensionless parameters as proxies for the light and the average quark mass as the physical point is approached.In the expressions of Φ 2 and Φ 4 , t 0 is the pion-and kaon-mass dependent flow observable; we use the notation t sym 0 to denote its value at the SU(3) f -symmetric point.We adopt the physical-point value √ 8t 0 = 0.4081(20)(37) fm from Ref. [71], obtained by equating the linear combination of pseudoscalar-meson decay constants to its physical value, set by the PDG values of the decay constants given below.Ref. [71] is an update of the work presented in [56] and includes a larger set of ensembles, including ensembles close to the physical point.We note that in Refs.[56,71] the absolute scale was determined assuming a slightly different definition of the physical point: the authors used the meson masses corrected for isospin-breaking effects as in [72], m π = 134.8(3)MeV and m K = 494.2(3)MeV.
Using the NLO χPT expressions, we have estimated the effect on f Kπ of these small shifts in the target pseudoscalar meson masses to be at the sub-permille level and therefore negligible for our present purposes.b.Strategy 2 Here we use f π -rescaling, which was already presented in our previous work [17], and express all dimensionful quantities in terms of the ratio f phys π /(af lat π ), where af lat π can be computed precisely on each ensemble.In this case, the intermediate scale t 0 is not needed and we use the following dimensionless proxies for the quark masses, As Φ 4 , the proxy y Kπ is approximately constant along our chiral trajectory.Since all relevant observables have been computed as part of this project, this method has the advantage of being fully self-consistent, and all correlations can be fully propagated.It will be our preferred strategy.
We use the following input to set the scale in our iso-symmetric theory [69,73], The quantity y Kπ is only used to correct for a small departure of the CLS ensembles from the physical value of this quantity, which we obtain using f K = 157.2(5)MeV [69,73].The latter, phenomenological value of f K implies a ratio f K /f π that is consistent with the latest lattice determinations [74][75][76].The impact of the uncertainty of f K on a win µ is small 2 , δa win µ 0.10×10 −10 , and occurs mainly through the strange contribution.In the isosymmetric theory, we take the phenomenological values of the triplet (m π , m K , f π ) as part of the definition of the target theory, and therefore only include the uncertainty from f K in our results.By contrast, in the final result including isospin-breaking effects, which we compare to a data-driven determination of a win µ , we include the experimental uncertainties of all quantities used as input.
The observables m π , m K , f π and f K , as well as t 0 /a 2 have been computed on all gauge ensembles and corrected for finite-size effects [77].Their values for all ensembles are listed in Table VII.

B. Fitting procedure
We now present our strategy to extrapolate the data to the physical point in our iso-symmetric setup.The ensembles used in this work have been generated such that the physical point is approached keeping 2 The sensitivity of a win µ to the value of fK can be derived from Table II.
approximately constant, where the two entries correspond respectively to Strategy 1 and 2.
To account for the small mistuning, only a linear correction in ∆X K = X phys K − X K is thus considered.To improve the fit quality, a dedicated calculation of the dependence of a win µ on X K has been performed, which is described in Appendix A. This analysis does not yet include all ensembles in the final result, and hence we decided to not apply this correction ensemble-byensemble prior to the global extrapolation to the physical point.Instead, we have used ∆X K to fix suitable priors on the fit parameter γ 0 in Eq. ( 26), which parametrizes the locally linear dependence on X K .The values of these priors are given in Appendix A.
To describe the light quark dependence beyond the linear term in (respectively for Strategy 1 and 2), we allow for different fit ansätze encoded in the function f ch (X π ).The precise choice of f ch is motivated on physical grounds and depends on the quark flavor.The specific forms will be discussed below.Since on-shell O(a)-improvement has been fully implemented, leading discretization artefacts are expected to scale as a 2 /t 0 up to logarithmic corrections [78,79].In the case of the vacuum polarization function, a further logarithmic correction proportional to a 2 log a was discovered in [80].Contrary to standard logarithmic corrections, it does not vanish as the coupling g 0 goes to zero due to correlators being integrated over very short distances.However, the intermediate window strongly suppresses the short-distance contribution, so that we do not expect this source of logarithmic enhancement to be relevant here.However, in the absence of further information on the relevant exponents of log a in full QCD [79], we still consider a possible logarithmic correction with unit exponent.Moreover, to check whether we are in the scaling regime, we consider higher order terms proportional to a3 .Finally, we also allow for a term ∝ X 2 a X π that describes pion-mass dependent discretization effects of order a 2 .
Thus, for each discretization of the vector correlator, the continuum and chiral extrapolation is done independently assuming the most general functional form where 'f' can be any flavor content and X a = a/ √ t 0 parametrizes the lattice spacing.Despite the availability of data from six lattice spacings and more than twenty ensembles, trying to fit all parameters is not possible.Thus each analysis is duplicated by switching on/off the parameters β 3 , δ and that control the continuum extrapolation.In addition, for each functional form f ch of the chiral dependence, different analyses are performed by imposing cuts in the pion mass (no cut, < 400 MeV, < 300 MeV) and/or in the lattice spacing.
Since several different fit ansätze can be equally well motivated, we apply the model averaging method presented in [81,82] where the Akaike Information Criterion (AIC) is used to weight different analyses and to estimate the systematic error associated with the fit ansatz (see also [20,83]).Thus, to each analysis (n) described above (defined by a specific choice of f ch , applying cuts in the pion mass or in the lattice spacing, and including or excluding terms proportional to β 3 , δ, ) we associate a weight w n given by where χ 2 is the minimum value of the chi-squared of the correlated fit, k is the number of fit parameters and n is the number of data points included in the fit 3 .The normalization factor N is such that the sum over all the analyses' weights are equal to one.Each analysis is again duplicated by either using the local-local or the local-conserved correlators.For those analyses, we use a flat weight.Finally, when cuts are performed, some fits may have very few degrees of freedom, and hence we exclude all analyses that contain fewer than three degrees of freedom.
The central value of an observable O is then obtained by a weighted average over all analyses and our estimate of the systematic error associated with the extrapolation to the physical point is given by The statistical error is obtained from the Jackknife procedure using the estimator defined by Eq. ( 28).

C. The continuum extrapolation at the SU(3) f -symmetric point
To reach sub-percent precision, a good control over the continuum limit is mandatory [79,80].As discussed below, it is one of the largest contributions to our total error budget.Thus, before presenting our final result at the physical point, we first demonstrate our ability to perform the continuum extrapolation.We have implemented three different checks: First, two discretizations of the vector correlator are used and the extrapolations to the physical point are done independently.Both discretizations are expected to agree within errors in the continuum limit.Physical observables computed using Wilson-clover quarks approach the continuum limit with a rate ∝ a 2 once the action and all currents are non-perturbatively O(a)-improved [53].
To check our ability to fully remove O(a) lattice artefacts in the action and the currents, two independent sets of improvement coefficients are used: both of them should lead to an a 2 scaling behavior but might differ by higher-order corrections.Finally, we have included six lattice spacings at the SU(3) f -symmetric point, all of them below 0.1 fm and down to 0.039 fm, to scrutinize the continuum extrapolation.In this section, we discuss those three issues, with a specific focus on the ensembles with SU(3) f symmetry.
Ensembles with six different lattice spacings in the range [0.039 : 0.099] fm are available for m π = m K ≈ 420 MeV.Since the pion masses do not match exactly, we first describe our procedure to interpolate our SU(3) f -symmetric ensembles to a single value of X π = X π , to be be able to focus solely on the continuum extrapolation.This reference point X π is chosen to minimize the quadratic sum of the shifts δX π = X π − X π .
We start by applying the finite-size effect correction discussed in the previous section to all ensembles.Then, a global fit over all the ensembles and simultaneously over both discretizations of the correlation function is performed using the functional form of Eq. ( 26) without any cut in the pion mass.Thus (γ 0 , γ 1 , γ 2 ) are fit parameters common to both discretizations, while the others are discretization-dependent.For the isovector contribution, we use the choice f ch (X π ) = 1/X π that leads to a reasonable χ 2 /d.o.f.= 1.1.The good χ 2 , and more importantly the good description of the light-quark mass dependence, ensures that the small interpolation to X π is safe and that we do not bias the result.In practice, we have checked explicitly that using different functional forms f ch to interpolate the data leads to changes that are small compared to the statistical error.Thus, for both choices of the improvement coefficients (set 1 and set 2), and for both discretizations LL and CL, the data from an SU(3) f -symmetric ensemble is corrected in the pseudoscalar masses to the reference SU(3) f -symmetric point at the same lattice spacing.The correction is obtained by taking the difference of Eq. ( 26) evaluated with the reference-point arguments (X a , X π , X K ) and the ensemble arguments (X a , X π , X K ), resulting in  where α = (LL), (CL) stands for the discretization.Note that X K = X π and X K = X π in view of the SU(3) f -symmetry.Throughout this procedure, correlations are preserved via the Jackknife analysis.
In a second step, we extrapolate both discretizations of the correlation function to a common continuum limit, using data at all six lattice spacings and assuming a polynomial in the lattice spacing, The two data sets obtained using the two different sets of improvement coefficients are fitted independently.The results are displayed in Fig. 2 for two cases: either applying f π -rescaling (left panel) or using t 0 to set the scale (right panel).For Set 1 of improvement coefficients, we observe a remarkably linear behavior over the whole range of lattice spacings, whether f πrescaling is applied or not.The second set of improvement coefficients (Set 2) leads to some visible curvature, but the continuum limit is perfectly compatible provided that lattice artefacts of order a 3 are included in the fit.We also tested the possibility of logarithmic corrections assuming the ansatz which is shown as the red symbol and red dashed curve in Fig. 2. The result is again compatible with the naive a 2 scaling, albeit with larger error.We conclude that logarithmic corrections are too small to be resolved in the data.We also remark that it is difficult to judge the quality of the continuum extrapolation based solely on the relative size of discretization effects between our coarsest and finest lattice spacing, as this measure strongly depends on the definition of the improvement coefficients.We tested the modification of the continuum extrapolation via X 2 a → (α s (1/X a )) ΓX 2 a as proposed in Refs.[79,84] for a win,I1 µ and a win,I0 µ ,c / in our preferred setup, using f π -rescaling and set 1 of improvement coefficients.The strong coupling constant α s has been obtained from the three-flavor Λ parameter of Ref. [85].Several choices of Γ in the range from 0.76 to 3 were tested.The curvature that is introduced by this modification, especially for larger values of Γ, would lead to larger values of a win µ in the continuum limit.However, such curvature is not supported by the data, as indicated by a deterioration of the fit quality when Γ is increased.Therefore, only small weights would be assigned to such fits in our model averaging procedure, where the modification has not been included.

D. Results for the isospin and flavor decompositions
Having studied the continuum limit at the SU(3) f -symmetric point, we are ready to present the result of the extrapolation to the physical point.The charm quark contribution is not included here and will be considered separately in Section V.
For the isovector or light quark contribution we use the same set of functional forms as in [17] The data shows some small curvature close to the physical pion mass.Thus, the variation f ch = 0 is excluded as it would significantly undershoot our ensemble at the physical pion mass (E250).We use Set 1 of improvement coefficients as our preferred choice and will use Set 2 only as a crosscheck.A typical extrapolation using f ch ( y) = 1/ y without any cut in the data is shown in the left panel of Fig. 3.We find that the specific functional form of f ch has much less impact on the extrapolation as compared to the inclusion of higher-order lattice artefacts.For the isoscalar and strange quark contributions, we restrict ourselves to functions that are not singular in the chiral limit: Again, the extrapolation using f ch ( y) = y log y with δ = 0 and without any cut in the data is shown in the right panel of Fig. 3.
Using the fit procedure described above, the AIC estimator defined in Eq. ( 28) leads to the following results for the isovector (I = 1) and the isoscalar contribution, charm excluded, where the first error is statistical and the second is the systematic error from the fit form used to extrapolate our data to the physical point.In Table II , we also provide the derivatives to translate our result to a different iso-symmetric scheme.We also note that both discretizations of the vector correlator yield perfectly compatible results.For the isovector contribution, and in units of 10 −10 , we obtain 186.14(0.87)stat (1.29) syst for the local-local discretization and 186.47(0.79)stat (0.79) syst for the local-conserved discretization, with a correlated difference of −0.33(0.72).For the isoscalar contribution, we find 47.39(0.24)stat (0.36) syst for the local-local discretization and 47.43(0.20)stat (0.19) syst for the local-conserved discretization, with a correlated difference of −0.04(0.10).
As an alternative to the fit weights given by Eq. ( 27), we have tried applying the weight factors used in Ref. [20]; see the footnote below Eq. (27).While a major change occurs in the subset of fits that dominate the weighted average, the results do not change significantly.
In particular, the central value of the isovector contribution changes by no more than half a standard deviation.
Finally, we have also performed an extrapolation to the physical point using the second set of improvement coefficients.Since our study at the SU(3) f -symmetric point shows curvature in the data, we exclude those continuum extrapolations that are only quadratic in the lattice  No rescaling f π rescaling FIG.4: Comparison of the isovector and isoscalar contributions (without the charm) using different variations (either using f π or t 0 to set the scale, and with both sets of improvement coefficients).The blue point is our final estimate obtained from the rescaling method with the set 1 of improvement coefficients.
spacing.The other variations are kept identical to those used for the first set.The results are slightly larger but compatible within one standard deviation.A comparison between the two strategies to set the scale and the two sets of improvement coefficients is shown in Fig. 4 for both the isovector and isoscalar contributions.
In order to facilitate comparisons with other lattice collaborations, we also present results for the light, strange and disconnected contributions separately.For the light and strange-quark connected contributions, we obtain For the disconnected contribution, the correlation function is very precise in the time range relevant for the intermediate window, and a simple sum over lattice points is used to evaluate Eq. ( 5).The data are corrected for finite-size effects using the method described in Section C. Since our ensembles follow a chiral trajectory at fixed bare average quark mass, we can consider a win,disc µ as being, to a good approximation, a function of the SU(3 Kπ )} (respectively for Strategy 1 and 2), with the additional constraint that the disconnected contribution vanishes quadratically in ∆ 2 for ∆ 2 → 0. We apply the following ansatz The ensembles close to the SU(3) f symmetric point (m π ≈ 350 MeV) are affected by significant FSE corrections and are not included in the fit.We obtain for the disconnected contribution and the extrapolation is shown in Fig. (5).The extrapolation using t 0 to set the scale shows less curvature close to the physical point.We use half the difference between the two extrapolations TABLE II: Derivatives of the window quantity a win µ (in units of 10 −10 ), for both the isovector and isoscalar contributions, as defined by Eq. (35).as our estimate for the systematic error.It is worth noting that the value for the intermediate window represents roughly 6% of the total contribution to a hvp,disc µ .As a crosscheck, we note that using Eqs.( 36), (37) and (39) we would obtain a win,I0 µ ,c / = (47.57±0.20 stat ±0.26 syst )×10 −10 , in good agreement with Eq. (34).

V. THE CHARM QUARK CONTRIBUTION
In our calculation, charm quarks are introduced in the valence sector only.A model estimate of the resulting quenching effect is provided in Appendix D. The method used to tune the mass of the charm quark has previously been described in Ref. [17] and has been applied to additional ensembles in this work.We only sketch the general strategy here, referring the reader to Ref. [17] for further details.For each gauge ensemble, the mass of the ground-state cs pseudoscalar meson is computed at four values of the charm-quark hopping parameter.Then the value of κ c is obtained by linearly interpolating the results in 1/κ c to the physical D s meson mass m Ds = 1968.35(0.07)MeV [69].We have checked that using either a quadratic fit or a linear fit in κ c leads to identical results at our level of precision.The results for all ensembles are listed in the second column of Table X.
The renormalization factor Ẑ(c) V of the local vector current has been computed nonperturbatively on each individual ensemble by imposing the vector Ward-identity using the same setup as in Ref. [58], but with a charm spectator quark.To propagate the error from the tuning of κ c , both Ẑ(c) V and a win,c µ are computed at three values of κ close to κ c .In the computation of correlation functions, the same stochastic noises are used to preserve the full statistical correlations.For both quantities, we observe a very linear behavior and a short interpolation to κ c is performed.The systematic error introduced by the tuning of κ c is propagated by computing the discrete derivatives of both observables with respect to κ c (second error quoted in Table X).This systematic error is considered as uncorrelated between different ensembles.
From ensembles generated with the same bare parameters but with different spatial extents (H105/N101 or H200/N202), it is clear that FSE are negligible in the charm-quark contribution.As in our previous work [17], the local-local discretization exhibits a long continuum extrapolation with discretization effects as large as 70% between our coarsest lattice spacing and the continuum limit, compared to only 12% for the local-conserved discretization.Thus, we discard the local-local discretization from our extrapolation to the physical point, which assumes the functional form Lattice artefacts are described by a polynomial in X a = a/ t sym 0 and a possible logarithmic term is included; recall that t sym 0 denotes the value of the flow observable at the SU(3) fsymmetric point.Only the set of proxies X π = φ 2 and X K = φ 4 is used.The light-quark dependence shows a very flat behavior, and a good χ 2 /d.o.f.= 0.9 is obtained without any cut in the pion mass.The corresponding extrapolation is shown on the right panel of Fig. 6.
Before quoting our final result, we provide strong evidence that our continuum extrapolation is under control by looking specifically at the SU(3) f -symmetric point where six lattice spacings are available.As for the isovector contribution, we use Eq. ( 40) to correct for the small pionmass mistuning at the SU(3) f -symmetric point.The data are interpolated to a single value of X * π using the same strategy as in Eq. (30).Those corrected points are finally extrapolated to the continuum limit using the ansatz (31).The result is shown in the left panel of Fig. 6 for the two sets of improvement coefficients of the vector current.Again, excellent agreement is observed between the two data sets.Even for the charm-quark contribution, we observe very little curvature when using the set 1 of improvement coefficients.
Having confirmed that our continuum extrapolation is under control, we quote our final result for the charm contribution obtained using the ansatz (40).Using Eq. ( 28), the AIC analysis described above leads to a win,c µ = (2.89± 0.03 stat ± 0.03 syst ± 0.13 scale ) × 10 −10 , where variations include cuts in the pion masses and in the lattice spacing, and fits where the parameters β 3 , β 4 and δ have been either switched on or off.

VI. ISOSPIN BREAKING EFFECTS
As discussed in the previous Sections III and IV A, our computations are performed in an isospin-symmetric setup, neglecting the effects due to the non-degeneracy of the up-and downquark masses and QED.At the percent and sub-percent level of precision it is, however, necessary to consider the impact of isospin-breaking effects.To estimate the latter, we have computed a win µ in QCD+QED on a subset of our isospin-symmetric ensembles using the technique of Monte Carlo reweighting [86][87][88][89][90] combined with a leading-order perturbative expansion of QCD+QED around isosymmetric QCD in terms of the electromagnetic coupling e 2 as well as the shifts in the bare quark masses ∆m u , ∆m d , ∆m s [90][91][92][93][94]. Consequently, we must evaluate additional diagrams that represent the perturbative quark mass shifts as well as the interaction between quarks and photons.We make use of non-compact lattice QED and regularize the manifest IR divergence with the QED L prescription [95], with the boundary conditions of the photon and QCD gauge fields chosen in accordance [93].We characterize the physical point of QCD+QED by the quantities m 2 π 0 , m 2 π 0 and the fine-structure constant α [91].The first three quantities are inspired by leading-order chiral perturbation theory including leading-order mass and electromagnetic isospin-breaking corrections [67], and correspond to proxies for the average light-quark mass, the strange-quark mass, and the lightquark mass splitting.As we consider leading-order effects only, the electromagnetic coupling does not renormalize [90], i.e. we may set e 2 = 4πα.The lattice scale is also affected by isospin breaking, which we however neglect at this stage.Making use of the isosymmetric scale [56], we match m 2 π 0 and m 2 K + + m 2 K 0 − m 2 π + in both theories on each ensemble and set π 0 to its experimental value.We have computed the leading-order QCD+QED quark-connected contribution to a win µ as well as the pseudoscalar meson masses m π 0 , m π + , m K 0 and m K + required for the hadronic renormalization scheme on the ensembles D450, N200, N451 and H102, neglecting quark-disconnected diagrams as well as isospin-breaking effects in sea-quark contributions.The considered quarkconnected diagrams are evaluated using stochastic U(1) quark sources with support on a single timeslice whereas the all-to-all photon propagator in Coulomb gauge is estimated stochastically by means of Z 2 photon sources.Covariant approximation averaging [96] in combination with the truncated solver method [97] is applied to reduce the stochastic noise.We treat the noise problem of the vector-vector correlation function at large time separations by means of a reconstruction based on a single exponential function.A more detailed description of the computation can be found in Refs.[91,92,98].The renormalization procedure of the local vector current in the QCD+QED computation is based on a comparison of the local-local and the conservedlocal discretizations of the vector-vector correlation function and hence differs from the purely isosymmetric QCD calculation [58] described in Section III B. We therefore determine the relative correction by isospin breaking in the QCD+QED setup.For f π -rescaling as introduced in Section IV A, isospin-breaking effects in the determination of f π are neglected.We observe that the size of the relative first-order corrections for a win µ is compatible on each ensemble and can in total be estimated as a (0.3 ± 0.1)% effect.

VII. FINAL RESULT AND DISCUSSION
We first quote our final result a win,iso µ in our iso-symmetric setup as defined in Section IV A. Using the isospin decomposition, and combining Eqs.(33), (34) where the first error is statistical, the second is the systematic error, and the last error of a win,iso µ is an estimate of the quenching effect of the charm quark derived in Appendix D. Overall, this uncertainty has a negligible effect on the systematic error estimate.The small bottom quark contribution has been neglected.For a hvp µ , this contribution has been computed in [99] and found to be negligible at the current level of precision.
As stressed in Section IV A, our definition of the physical point in our iso-symmetric setup is scheme dependent.To facilitate the comparison with other lattice collaborations, the derivatives with respect to the quantities used to define our iso-symmetric scheme are provided in Table II.They can be used to translate from one prescription to another a posteriori.
One of the main challenges for lattice calculations of both a hvp µ and the window observable is the continuum extrapolation of the light quark contribution, which dominates the results by far.To address this specific point, we have used six lattice spacings in the range [0.039,0.0993]fm in our calculation, along with two different discretizations of the vector current (see the discussion in Section IV C).Although this work contains many ensembles away from the physical pion mass, we observe only a mild dependence on the proxy used for the light-quark mass.This observation is corroborated by the fact that, in the model averaging analysis, most of the spread comes from fits that differ in the description of lattice artefacts rather than on the functional form f ch that describes the light-quark mass dependence.In Fig. 7, we compare our results in the isosymmetric theory with other lattice calculations.Our estimate for a win,iso µ agrees well with that of the BMW collaboration who quote a win,iso µ = 236.3(1.4) × 10 −10 using the staggered quark formulation [20].However, our result is about 2.3σ above the published value by the RBC/UKQCD collaboration, a win,iso µ = 232.0(1.5)× 10 −10 , obtained using domain wall fermions [13].It is also 1.7σ above the recent estimate quoted by ETMC, based on the twisted-mass formalism [22], which reads a win,iso µ = 231.0(2.8)×10−10 .The difference with the latter two calculations can be traced to the light-quark contribution a win,ud µ , which is shown in the second panel from the right.In this context, it is interesting to note that, apart from BMW, two independent calculations using staggered quarks (albeit with a different action as compared to the BMW collaboration) have quoted results for a win,ud µ [18,21,24] that are in good agreement with our estimate, as can be seen in Fig. 7.The middle panel of the figure shows that our estimate for the strange quark contribution is slighly higher compared to other groups, but due to the relative smallness of a win,s µ this cannot account for the difference between our result for a win,iso µ and Refs.[22] and [13].Good agreement with the BMW, ETMC and RBC/UKQCD collaborations is found for both the charm and quark-disconnected contributions.
If one accepts that most lattice estimates for the light-quark connected contribution a win,ud µ have stabilized around ≈ 207 × 10 −10 , one may search for an explanation why the results by RBC/UKQCD [13] and ETMC [22] are smaller by about 2%.This is particularly important since a win,ud µ contributes about 87% to the entire intermediate window observable.One possibility is that the extrapolations to the physical point in Refs.[13] and [22] are both quite long.For instance, the minimum pion mass among the set of ensembles used by ETMC is only about 220 MeV, while the result by RBC/UKQCD has been obtained from two lattice spacings, i.e. 0.084 fm and 0.114 fm.Further studies using additional ensembles at smaller pion mass and lattice spacings are highly desirable to clarify this important issue.
In order to compare our result with phenomenological determinations of the intermediate window observable, we must correct for the effects of isospin-breaking.Our calculation of isospinbreaking corrections, described in Section VI, has been performed on a subset of our ensembles and is, at this stage, lacking a systematic assessment of discretization and finite-volume errors.Furthermore, only quark-connected diagrams have been considered so far.To account for this source of uncertainty, we double the error and thereby apply a relative isospin-breaking correction of (0.3 ± 0.2)% to a win,iso µ , which amounts to a shift of +(0.70 ± 0.47) × 10 −10 .Thus, our final result including isospin-breaking corrections is a win µ = (237.30± 0.79 stat ± 1.13 syst ± 0.05 Adding all errors in quadrature yields 237.30(1.46)× 10 −10 which corresponds to a precision of 0.6%.A comparison with other lattice calculations is shown in Fig. 8. Since corrections due to isospin breaking are small, the same features are observed as in the isosymmetric theory: including isospin-breaking corrections with the estimates by ETMC [22], BMW [20] and RBC/UKQCD [13].The estimate based on the data-driven method of Ref. [48] is shown in red.
while our result agrees well with the published estimate from BMW [20], it is larger than the values quoted by ETMC [22] and RBC/UKQCD [13].Our result lies 3.9σ above the recent evaluation using the data-driven method [48], which yields a win µ = 229.4(1.4) × 10 −10 and is shown in red in Fig. 8.Our result for a win µ is also consistent with the observation that the central value of our 2019 result for the complete hadronic vacuum polarization contribution [17] lies higher than the phenomenology estimate, albeit with much larger uncertainties.In Ref. [62] we observed a similar, but statistically much more significant enhancement in the hadronic running of the electromagnetic coupling, ∆α had (−Q 2 ) relative to the data-driven evaluation, especially for Q 2 < ∼ 3 GeV 2 .As pointed out at the end of Section II, the relative contributions from the three intervals of center-of-mass energy separated by √ s = 600 MeV and √ s = 900 MeV are similar for a win µ and ∆α had (−1GeV 2 ), even though the respective weight functions in the time-momentum representation are rather different.The fact that the lattice determination is larger by more than three percent for both quantities, in each case with a combined error of less than one percent, suggests that a genuine difference exists at the level of the underlying spectral function, R(s)/(12π 2 ), between lattice QCD and phenomenology.
If one were to subtract the data-driven evaluation of a win µ from the White Paper estimate [3] and replace it by our result in Eq. ( 45), the tension between the SM prediction for a µ and experiment would be reduced to 2.9σ.This observation illustrates the relevance of the window observable for precision tests of the SM.Our findings also strengthen the evidence supporting a tension between data-driven and lattice determinations of a hvp µ .In our future work we will extend the calculation to other windows and focus on the determination of the full hadronic vacuum polarization contribution, a hvp µ .
programme, AMX-18-ACE-005 and from the French National Research Agency under the contract ANR-20-CE31-0016.We are grateful to our colleagues in the CLS initiative for sharing ensembles.
Appendix A: Mistuning of the chiral trajectory The ensembles used in our work have been generated with a constant bare average sea quark mass which differs from a constant renormalized mass by O(a) cutoff effects.When the sum of the renormalized quark masses is kept constant, the dimensionless parameters φ 4 and y Kπ , which have been introduced in Section IV A to define the chiral trajectories towards the physical point, are constant to leading order in chiral perturbation theory (χPT).Therefore, φ 4 and y Kπ cannot be constant across our set of ensembles due to cutoff effects and higher-order effects from χPT.
We have to correct for the sources of mistuning of our ensembles with respect to the chiral trajectories of strategies 1 and 2. This can be done by parameterizing the dependence of our observables on X K ∈ {y Kπ , φ 4 } in the combined chiral-continuum extrapolation.However, since the pion and kaon masses are not varied independently within our set of ensembles, the dependence on ∆X K = X phys K − X K cannot be resolved reliably in our fits.A different strategy has to be employed to stabilize our extrapolation to the physical point.
Explicit corrections of the mistuning prior to the chiral extrapolation have been used in [56] to approach the physical point at constant φ 4 = φ phys 4 .These corrections are based on small shifts defined from the first order Taylor expansion of the quark mass dependence of lattice observables.The expectation value of a shifted observable is given by with the N f = 3 sea quark mass shifts ∆m q,i .Within this appendix, we work with observables and expectation values that are defined after integration over the fermion fields, i.e. the expectation values are taken with respect to the gauge configurations.The total derivative of an observable with respect to the quark masses is decomposed via The partial derivative of an observable with respect to a quark mass of flavor i captures the effect of shifts of valence quark masses.The second and third terms that contain the derivative of the action S with respect to the quark masses account for sea quark effects.The chain rule is used to compute the derivatives of derived observables.The chain rule relating the derivatives with respect to the quark masses to those with respect to the variables X j = X π , X K can be written ∀ n = (n 1 , n 1 , n 3 ), the condition n 1 = n 2 being imposed to remain in the isosymmetric theory.
In particular, if the direction of the vector n in the space of quark masses is chosen such that ∆ π ( n) vanishes, the following expression [71] for the derivative of an observable with respect to In [56] the shifts n i have been chosen to be degenerate for all three sea quarks.In [71] the same approach is taken at the SU(3) f -symmetric point and n = (0, 0, 1) is used when am q,l = am q,s .
To stabilize the predictions for the derivatives, they are modeled as functions of lattice spacing and quark mass.
To improve the reliability of our chiral extrapolation, we have determined the derivatives of a win,ud µ and a win,s µ with respect to light and strange quark masses on a large subset of the ensembles in Table I.Whereas the computation of the first term in Eq. (A2) shows a good signal for the vector-vector correlation function, the second and third term carry significant uncertainties.In the case of f π -rescaling, a non-negligible statistical error that has its origin in df π /dm q,i enters the derivative of a win µ .Our computation does not yet cover all ensembles in this work and has significant uncertainties on some of the included ensembles.Moreover, we have not computed the mass-derivative of a win,disc µ that enters a win,I0 µ .Therefore, we have decided not to correct our observables prior to the global extrapolation but to determine the coefficient γ 0 in Eq. ( 26) instead.We do not aim for a precise determination here but focus instead on the determination of a sufficiently narrow prior width, in order to stabilize the chiral-continuum extrapolation.
We compute the derivatives with respect to X K as specified in Eq. (A4) with the shift vector n chosen such that ∆ π ( n) vanishes ensemble by ensemble, i.e. the shift is taken in a direction in the quark mass plane where X π remains constant.The derivatives are therefore sensitive to shifts in the kaon mass.A residual shift of X a is present at the permille level.
We collect our results for the derivatives with respect to φ 4 and y Kπ in Table III.Throughout this appendix, we use units of 10 −10 for a win µ , as well as for coefficient γ 0 .The results are based on the local-local discretization of the correlation functions and the improvement coefficients and renormalization constants of set 1.As can be seen, the derivative of the isovector contribution to the window observable vanishes within error on most of the ensembles.This is expected from the order-of-magnitude estimate in Eq. (B33).No clear trend regarding a dependence on X π , X K or X a can be resolved.We show the derivative of a win,I1 µ with respect to X π in the upper panels of Fig. 9.For the corresponding priors for the chiral-continuum extrapolation we choose The derivative of the strange-connected contribution of the window observable with respect to X K is negative and can be determined to good precision.Our results are shown in the lower panels of Fig. 9.We choose our priors such that their width encompasses the spread of the data.For the strange-connected and the isoscalar contribution, we choose γ win,s,y Kπ These values are compatible with the estimate in Eq. (B26).Discretization effects in the data may be inspected by comparing the derivatives based on the two sets of improvement coefficients.Such effects are largest for the two ensembles at β = 3.34, but are still smaller than the spread in the data and therefore not significant with respect to our prior widths.In our global extrapolations, we use a single set of priors irrespective of the improvement procedure.

Sensitivity of the window quantity
In [49], a semi-realistic model for the R-ratio was used for the sake of comparisons with lattice data generated in the (u, d, s) quark sector with exact isospin symmetry.In particular, the model does not include the charm contribution, nor final states containing a photon, such as π 0 γ.It leads to the following values for the window observables and their sum, the full a hvp µ , Given the omission of the aforementioned channels, these values are quite realistic. 4Here we only use the model to provide the partition of the quantities above into three commonly used intervals of √ s, in order to illustrate what the relative sensitivities of these quantities are to different energy intervals.These percentage contributions are given in Table IV, along with the corresponding figures for the subtracted vacuum polarization, The model yields for this quantity the value 385.5 × 10 −4 at Q 2 = 1 GeV 2 .We expect the fractions in the table to be reliable with an uncertainty at the five to seven percent level.and the partial quantities (a hvp µ ) SD,ID,LD , as well as the subtracted vacuum polarization at scale Q 2 = 1 GeV 2 , according to the R-ratio model given in [49].Note that this model includes neither the charm nor final states containing a photon, such as π 0 γ.
The model value for the intermediate window is best compared to the sum of Eqs.(33,34).The difference is (1.8 ± 1.4) × 10 −10 , which represents agreement at the 1.3σ level.The main reason the R-ratio model agrees better with the lattice result than a state-of-the-art analysis [48] is that the model does not account for the strong suppression of the experimentally measured R-ratio in the region 1.0 < √ s/GeV < 1.5 relative to the parton-model prediction.This observation suggests a possible scenario where the higher lattice value of a win µ as compared to its data-driven evaluation is explained by a too pronounced dip of the R-ratio just above the φ meson mass.In such a scenario, the relative deviation between the central values of a hvp µ obtained on the lattice and using e + e − data would be smaller than for a win µ by a factor of about 1.5, given the entries in Table IV.Indeed, it has been shown [50] that the central values of the BMW collaboration [20] cannot be explained by a modification of the experimental R(s) ratio below s = 1 GeV 2 alone.

Model estimate of
In [62], we have used two closely related R-ratio models for the strangeness correlator and the light-quark contribution to the isoscalar correlator, with m ω = 0.78265 GeV, m φ = 1.01946GeV and [100] The threshold values s 0 and s 1 have been adjusted to reproduce the corresponding lattice results for a hvp µ .The model R-ratios of Eqs.(B8-B9) were used [62] in the linear combination (18R I=0 − 9R s ) in order to model the SU(3) f breaking contribution Π 08 , which enters the running of the electroweak mixing angle.Our model for this linear combination also obeys an exact sum rule, ∞ 0 ds (18R I=0 − 9R s ) = 0, within the statistical uncertainties.We now evaluate the window quantity for the models of Eqs.(B8-B9).For the strangeness contribution, we have Given the modelling uncertainties, these values are in excellent agreement with the lattice results presented in the main part of the text, respectively Eqs. (37) and (34).We also record some useful values of the kernel, In the following, we evaluate the strange-quark mass dependence of a win,s µ , based on the idea that the parameters A φ , m φ and s 1 only depend on the mass of the valence (strange) quark.This general assumption is reflected in Eqs.(B21, B22, B24) below.
It was noted a long time ago [101] that the electronic decay widths of vector mesons, normalized by the relevant charge factor, is only very weakly dependent on their mass: This suggests that, unlike in QED, (A V • m V ) depends less strongly on m V than A V itself for QCD vector mesons.Therefore it is best to estimate the derivative of interest as follows, We estimate the following derivatives by taking a finite difference between the ω and the φ meson properties, and To the statistical error from the electronic widths of the ω and φ mesons, we have added a modelling error of 15%.Using t 0 , the value above translates into ∂a win,s µ ∂φ 4 φ 2 (−10.9 ± 0.7 stat ± 1.6 model ) × 10 −10 , (B27) which can directly be compared to the values from lattice QCD listed in Table III.The agreement is excellent.In Eq. (B9), we have written the perturbative contribution above the threshold s 1 in the massless limit.We now verify that the mass dependence of the perturbative contribution is negligible for fixed s 1 .The leading mass-dependent perturbative contribution to the R-ratio well above threshold is (see e.g.[102], Eqs.11 and 12) From here we have estimated ∂ (B29) Here, the dependence on s 1 only contributes 18% of the total.We have again assigned a 15% modelling uncertainty to the prediction.Since we expect valence-quark effects to dominate, the prediction (B29) can also be applied to the full isoscalar prediction.
The influence of the strange quark mass on the isovector channel is a pure sea quark effect, and is as such harder to estimate.Based on the OZI rule, one would also expect a smaller relative sensitivity than in the strangeness channel addressed in the previous subsection.
One effect of the presence of strange quarks on the isovector channel is that kaon loops can contribute.No isovector vector resonances with a strong coupling to KK are known, therefore we attempt to use scalar QED (sQED) to evaluate the effect of the kaon loops.Note that at the SU(3) f symmetric point, the sum of the K0 K 0 and K + K − contributions to the isovector channel amounts to half as much as that of the pions.We find, integrating in s from threshold up to 4 GeV (B31) A further, more indirect effect of two-kaon intermediate states is that they can affect the properties of the ρ meson.On general grounds, one expects the two-kaon states to reduce the ρ mass, since energy levels repel each other.However, for the window quantity it so happens that sf win (s) has a maximum practically at the ρ mass, therefore the derivative of this function is extremely small, 2 f win (s) d ds (sf win (s)) The effect of a shift in the ρ meson mass is therefore heavily suppressed. 5Reasonable estimates of the order-of-magnitude of the derivative µ which is smaller than the sQED estimate.These estimates are based on the observation that the ratio m ρ /f π is about 5% higher at a pion mass of 311 MeV in the N f = 2 QCD calculation [103] than if one interpolates the corresponding N f = 2 + 1 QCD results [17,104] to the same pion mass, though a caveat is that neither result is continuum-extrapolated.The effect of the kaon intermediate states on the ππ line-shape is even harder to estimate, but we note that even in N f = 2 QCD calculations [103], i.e. in the absence of kaons, the obtained g ρππ coupling is consistent with N f = 2 + 1 QCD calculations [17,104] carried out at comparable pion masses.In summary, we use the sQED evaluation of Eq. (B31) to provide the order-of-magnitude estimate We note that the statistical precision of our lattice-QCD results for this derivative in Table III is not sufficient to resolve the small effect estimated here.
Appendix C: Finite-volume correction Corrections for finite-size effects (FSE) have been estimated using a similar strategy to the one presented in our previous publication on the hadronic contributions to the muon g − 2 [17].The main difference lies in the treatment of small Euclidean times, where we have replaced NLO χPT by the Hansen-Patella method as described below.We have also investigated finite-size corrections in χPT at NNLO [20,105].Overall, we found it to be comparable in size to the values found in Tables V-VI, the level of agreement improving for increasing volumes and decreasing pion masses.Given that the NNLO χPT correction term is in many cases not small compared to the NLO term, we refrain from using χPT to compute finite-size effects in our analysis of a win µ (see [24] for a more detailed discussion of the issue).

The Hansen-Patella method
In [106,107], finite-size effects for the hadronic contribution to the muon (g −2) are expressed in terms of the forward Compton amplitude of the pion as an expansion in exp (−| n|m π L) for | n| 2 = 1, 2, 3, 6, . . . .Here, n k schematically represents the number of times the pion propagates around the k th spatial direction of the lattice.Corrections that start at order exp (−n eff m π L) with n eff = 2 + √ 3 ≈ 1.93 are neglected: they appear when at least two pions propagate around the torus.The results for the first three leading contributions (| n| 2 ≤ 3) can thus be used consistently to correct the lattice data on each timeslice separately.We decided to use the size of the | n| 2 = 3 term, i.e. the last one that is parametrically larger than the neglected n eff ≈ 1.93 contribution, as an estimate of the inherent systematic error.
In this work we follow the method presented in [107], where the forward Compton amplitude is approximated by the pion pole term, which is determined by the electromagnetic form factor of the pion in the space-like region.Since the form-factor is only used to evaluate the small finite-volume correction, a simple but realistic model is sufficient.Here we use a monopole parametrization obtained from N f = 2 lattice QCD simulations [108], The statistical error on the finite-size correction is obtained by propagating the jackknife error on the pion and monopole masses.The results obtained using this method are summarized in the third and fourth columns of Table V and Table VI.

The Meyer-Lellouch-Lüscher formalism with Gounaris-Sakurai parametrization
As an alternative, we also consider the Meyer-Lellouch-Lüscher (MLL) formalism.The isovector correlator in both finite and infinite volume is written in terms of spectral decompositions where F π (ω) is the time-like pion form factor.Following the Lüscher formalism, the discrete energy levels E i = 2 m 2 π + k 2 i in finite volume are obtained by solving the equation where φ(q) is a known function [109,110], n a strictly positive integer and δ 1 is the scattering phase shift in the isospin I = 1, p-wave channel.Strictly speaking, this relation holds exactly only below the four-particle threshold that starts at 4m π .This is only a restriction at light pion mass where many states are needed to saturate the spectral decomposition in finite volume.We will see below how to circumvent this difficulty.In [111], the overlap factors A i that enter the spectral decomposition in finite volume were shown to be related to the form factor in infinite volume through the relation The time-like pion form factor has been computed on a subset of our lattice simulations [17,104].Since the form factor is only needed to estimate the small finite-volume correction, an approximate model can be used.Here, we assume a Gounaris-Sakurai (GS) parametrization that contains two parameters: the g ρππ coupling and the vector meson mass m ρ [112].A given choice of those parameters allows us to compute both the finite-volume and infinite-volume correlation function in the isovector channel at large Euclidean times using Eq.(C3).The difference G I=1 (t, ∞) − G I=1 (t, L), when inserted into Eq.( 5), yields our estimate of the FSE.
In practice, the GS parameters are obtained from a fit to the isovector correlation function G I=1 (t, L) at large Euclidean times, using Eqs.(C3), (C4) and (C5).Statistical errors on the GS parameters can easily be propagated using the Jackknife procedure.Since this method is expected to give a good description only up to the inelastic threshold, Eq. (C4) being formally valid below 4m π , we opt to use the MLL formalism only above a certain cut in Euclidean time, given by t * = (m π L/4) 2 /m π .Below the cut, we always use the HP method described above.Above the cut, the lightest few finite-volume states in the spectral decomposition saturate the integrand.The results using the MLL formalism are summarized in the fifth column of Table V and Table VI.

Corrections applied to lattice data
In Tables V and VI we summarize the FSE correction applied to the raw lattice data.We find that finite-size corrections computed using either the HP or the MLL method for (t > t ) show good agreement within their respective uncertainties.Our final estimates, shown in the rightmost column, are obtained by adding the result from the HP method at short times (t < t ) to that of the MLL method above t and the kaon loop contribution.The latter has been computed in χPT at NLO (see for instance [113]) on ensembles without SU(3) flavor symmetry.At the SU(3) symmetric point, the kaon loop contribution has been accounted for by scaling the HP and MLL corrections by a factor of 3/2.We have included the scale factor in the respective entries in Tables V and VI.
The uncertainty quoted in the rightmost column is given by the statistical error computed as described in the two previous sections.It includes the statistical error on the GS parameters and on the monopole mass that appears in the parametrization of the form factor in Eq. (C1).The systematic error on the HP contribution is estimated as described in Section C 1.
For our final estimates of finite-volume corrections, we adopt a more conservative approach regarding the overall uncertainty.As in our earlier paper [10], we base our uncertainty estimate on the comparison to the NLO χPT correction, which leads us to assign an error of 25% of the estimated correction for each ensemble, which replaces the uncertainties quoted in the last column of Tables V and VI.For example, the finite-size correction applied to a win,I1 µ in the case of ensemble J303 with f π -rescaling is (1.62 ± 0.405) × 10 −10 .TABLE V: Finite-size effects in the isovector channel with f π -rescaling, in units of 10 −10 , for our ensembles described in Table I.The correction obtained using the HP method is given in the third and fourth columns.The MLL estimate in the long-distance region is listed in the fifth column.The contribution of the kaon is given in column six, where dashes for ensembles at the SU (3)  The gauge configurations used in this work contain the dynamical effects of up, down and strange quarks.As for the charm quarks, we have only taken into account the connected valence contributions.In this appendix, we estimate the systematic error from the missing effect of charm sea-quark contributions.The question we are after can be formulated as, "What is the charm-quark effect on the R-ratio in a world in which the charm quark is electrically neutral?".
As in [62], we adopt a phenomenological approach.There, we evaluated the perturbative prediction for the charm sea quark effect and found it to be small for the running of the electromagnetic coupling from Q 2 = 1 GeV 2 to 5 GeV 2 .Alternatively, we considered D-meson pair creation in the electromagnetic-current correlator of the (u, d, s) quark sector.The contribution of the D + D − channel to the R-ratio reads and similar expressions hold for the D 0 D0 and D + s D − s channels.Since the form factor F D + is not known precisely and our goal is only to estimate the order of magnitude of the effect, we will approximate it by its value at s = 0, which amounts to treating D-mesons in the scalar QED framework and replacing their form factors by the relevant electromagnetic charges: Note that up-, down-, or strange-quarks play the role of the valence quarks giving the mesons their respective charges.The corresponding contributions to a hvp where we have inserted the a hvp µ = 720.0value from Ref. [17].The charm sea-quark contributions are thus negligible at the current level of precision.
We notice that ∆ c-sea a hvp µ /a hvp µ is much smaller than the effects found in the HVP contributions to the QED running coupling, namely ∼ 0.4% [62].We interpret the difference as follows: the typical scale in a hvp µ is given by the muon mass, which is well separated from the D-meson masses.Therefore the D-meson effects are strongly suppressed.In comparison, the running coupling was investigated at the GeV scale and the suppression is less strong.
In the intermediate window, the charm sea-quarks are even more suppressed, as seen in the tiny value of ∆ c-sea a win µ /a win µ .This results from the following fact: creating D-meson pairs requires a center-of-mass energy of ∼ 4 GeV, corresponding to t ∼ 0.05 fm, which is much smaller than the lower edge of the intermediate window, t 0 = 0.4 fm.Therefore, the D-meson pair creation contributes mostly to the short-distance window (a hvp µ ) SD .In fact, the effect in the intermediate window ∆ c-sea a win µ amounts to at most 5% of the total ∆ c-sea a hvp µ .Charm sea quarks lead not only to on-shell D mesons in the R(s) ratio, but also to virtual effects below the threshold for charm production.This is seen explicitly in the perturbative calculation [115], where the two effects are of the same order.At present, we do not have a means to estimate these virtual effects on the quantity a win µ , in which they are less kinematically suppressed.Therefore, we will conservatively amplify the uncertainty that we assign to the neglect of sea charm quarks by a factor of three relative to the prediction of Eq.D6.This estimate also generously covers the effect on a win µ which follows from adopting the perturbative charm-loop effect on R(s) down to s = 1 . . .1.5 GeV 2 .Thus, rounding the uncertainty to one significant digit, we quote ∆ c-sea a win µ = 0.05 × 10 −10 (D7) as the uncertainty on a win µ due to the quenching of the charm in the final result Eq. ( 44) for the isosymmetric theory.

FIG. 2 :
FIG.2: Continuum extrapolation for the isovector quark contribution at the SU(3) f -symmetric point.Left: using f π -rescaling.Right: with t 0 to set the scale.The blue and green points correspond to the two different sets of improvement coefficients (see Section III).For clarity, the extrapolated results have been shifted to the left.

85 FIG. 3 :
FIG.3: Left: one typical extrapolation of the isovector contribution using f ch ( y) = 1/ y.The data corresponds to the local-conserved discretization of the correlator using the set 1 of improvement coefficients.Error bands are the results from the fit for each of the six lattice spacings.The black line is the chiral extrapolation in the continuum limit.The black point is the result at the physical point.Right: same for the isoscalar contribution but using f ch ( y) = 0.

Kπ ) 2 a 70 FIG. 5 :
FIG.5: Extrapolation to the physical point for the quark-disconnected contribution using Eq.(38).The vertical dashed line represents the physical point in our iso-symmetric QCD setup.The black point is the result of the extrapolation, and the grey band represents the extrapolation to the continuum limit with X K = X K .Points with dashed error bars are not included in the fit.

a 2 [ 85 FIG. 6 :
FIG. 6: Left panel: study of the continuum extrapolation of the charm quark contribution to a win µ at the SU(3) f -symmetric point using the local-conserved discretization of the correlation function.The black and green points are obtained using two independent sets of improvement coefficients, as explained in Section III B. Right panel: Example of a typical extrapolation to the physical point of the charm-quark contribution.The error from the scale setting, which is highly correlated between ensembles, is not shown.The plain lines are obtained from the fit function (40) without any cut in the pion mass.

FIG. 7 :
FIG.7: Comparison of our results (in units of 10 −10 ) with other lattice calculations[13,18,[20][21][22][23][24] in isosymmetric QCD.The four panels on the left show compilations of the individual quark-disconnected, charm, strange and light quark contributions.The total result for a win µ in the isosymmetric case is shown in the rightmost panel.Our results are represented by green circles and vertical bands.

FIG. 8 :
FIG. 8: Comparison of our result for a win µ

85 FIG. 9 :
FIG.9: Derivatives of the isovector and the strange-connected contributions to the window observable with respect to X π .The gray areas illustrate the priors that are used in the global extrapolation.

∂m 2 K 2 K
a win,s,pert µ ≈ 0.5 × 10 −10 GeV −2 .Since this contribution to ∂ ∂m a win,s µ is about one sixth the statistical uncertainty from the vector meson electronic decay widths, we neglect the perturbative mass dependence of a win,s µ .For future reference, we evaluate in the same way as in Eq. (B26) the derivative of a hvp,s 129 ± 6 stat ± 19 model ) × 10 −10 GeV −2 .
hvp (s) R D 0 D 0 + R D + D − + R D +where m µ is the muon mass and the analytic form of K(s) can be found e.g.in[114], section 4.1.Similarly, the counterpart for the intermediate window reads∆ c-sea a win µ = ∞ 0 ds f win (s) R D 0 D 0 + R D + D − + R D + s D − s (s) .(D4)where f win (s) is defined in Eq. (B2).For the D-meson masses, we use the values provided by the Particle Data Group 2020[100].Our results are ∆ c-sea a hvp µ symmetric point indicate that this contribution is contained in the HP and MLL estimates.Our final estimate is given in the last column.Only statistical errors are shown.We assign an uncertainty of 25% of the FSE on each ensemble (see text).

TABLE VI :
Same as Table V using t 0 to set the scale.

TABLE X :
Charm hopping parameter κ c , renormalization factor of the local vector current Z for the two discretizations of the correlator: local-local (LL) and local-conserved (LC).The first error is statistical and the second is the systematic error arising from the tuning of the charm quark hopping parameter as explained in the main text.