Hadronic vacuum polarization: comparing lattice QCD and data-driven results in systematically improvable ways

The precision with which hadronic vacuum polarization (HVP) is obtained determines how accurately important observables, such as the muon anomalous magnetic moment, a_\mu, or the low-energy running of the electromagnetic coupling, \alpha, are predicted. The two most precise approaches for determining HVP are: dispersive relations combined with e+e- to hadrons cross-section data, and lattice QCD. However, the results obtained in these two approaches display significant tensions, whose origins are not understood. Here we present a framework that sheds light on this issue and, if the two approaches can be reconciled, allows them to be combined. Via this framework, we test the hypothesis that the tensions can be explained by modifying the R-ratio in different intervals of center-of-mass energy sqrt(s). As ingredients, we consider observables that have been precisely determined in both approaches. These are the leading hadronic contributions to a_\mu, to the so-called intermediate window observable and to the running of \alpha between spacelike virtualities 1GeV^2 and 10GeV^2 (for which only a preliminary lattice result exists). Our tests take into account all uncertainties and correlations, as well as uncertainties on uncertainties in the lattice results. Among our findings, the most striking is that results obtained in the two approaches can be made to agree for all three observables by modifying the \rho peak in the experimental spectrum. In particular, we find that this requires a common ~5\% increase in the contributions of the peak to each of the three observables. This finding is robust against the presence or absence of one of the constraining observables. However, such an increase is much larger than the uncertainties on the measured R-ratio. We also discuss a variety of generalizations of the methods used here, as well as the limits in the information that can be extracted...


I. INTRODUCTION
A virtual, propagating photon polarizes the vacuum into quarks and gluons.Known as hadronic vacuum polarization (HVP), this effect is important when processes involving electromagnetism at the quantum level are studied with high precision.For a photon of small virtuality, predicting this effect requires being able to describe the strong interaction in its nonperturbative regime.At present, there are two approaches for making precise predictions of this polarization.These are based, on the one hand, on large-scale, numerical simulations in lattice quantum chromodynamics (QCD); and on the other, on the exploitation of e + e − → hadrons data.
HVP has been recently the center of much attention, because of its importance in the standard model prediction of the anomalous magnetic moment of the muon, a µ .This quantity is currently being measured at the Fermi National Accelerator Laboratory (FNAL) [1] and was previously measured at the Brookhaven National Laboratory (BNL) [2].While the data-driven determination of the leading-order (LO) HVP contribution to a µ [3][4][5], a LO-HVP µ , yields a standard model prediction 4.2σ below the measurement of the total a µ [1], the most precise lattice calculation, of this contribution [6], reduces this difference to 1.5σ, with comparable uncertainties [7].As discussed below, and in Ref. [6], this reduction in the tension between prediction and experiment comes at the expense of being 2.0σ above the data-driven prediction.Moreover, for the so-called intermediate-window contribution to a LO-HVP µ [8], this excess rises to 3.8σ ( [6] and below), and even 4σ for a weighted average of independent lattice determinations of this quantity [6,[9][10][11][12] (see Sec. IV A).
The above discussion does not take into account the measurement of the e + e − → π + π − (γ) cross section, from threshold to 1.2 GeV, reported by the CMD-3 collaboration in the preprint [13].This is because an explanation has not yet been found for the fact that the cross section which they obtain is significantly larger than all previous, modern measurements and, in particular, the one published by the same collaboration in Ref. [14].
Another place where HVP plays an important role is in the scale-dependence of the electromagnetic coupling, α.The fact that the lattice predicts a larger value for a LO-HVP µ has an impact on the value of α at the scale of arXiv:2308.04221v1[hep-ph] 8 Aug 2023 the Z-boson mass, M Z .This has been investigated in Refs.[15][16][17][18].The overall conclusion is that the lattice excess in a µ , alone, does not imply a change in α(M2 Z ) large enough to have an impact on precision electroweak tests.This is confirmed by direct lattice calculations of the running of α [6,19,20].In particular, Ref. [6] suggests that the faster running of α observed in lattice calculations is concentrated for spacelike values of momentum scale below a few GeV 2 .
Beyond signaling tensions between the lattice and data-driven approaches for specific physical observables, the discrepancies discussed above also contain information about the agreement or disagreement between the primary quantities used to compute these HVP related observables in each of the two approaches.On the lattice this primary observable is the quark electromagnetic-current, two-point function with vanishing three-momentum, computed as a function of Euclidean time t.In the data-driven approach, it is instead the cross section for e + e − annihilation into hadrons, measured as a function of center-of-mass (c.o.m.) energy √ s, normalized by the tree-level cross section for e + e − → µ + µ − in the massless limit, i.e. the R-ratio R(s). 1o help pinpoint the possible sources of the disagreement between the two approaches, and possibly correct them, one needs to be able to make sharp statements about the agreement or disagreement between the lattice correlation function and the experimental R-ratio in different regions of t and √ s.Because the current correlator is proportional to the Laplace transform of √ sR(s) [21], it is straightforward to determine that correlator, and any observable that can be obtained from it, once R(s) is measured.Thus, it is relatively simple to identify regions of t in which the lattice and data-driven approaches disagree, possibly pointing to a problem with lattice computations at those length scales.
Determining the R-ratio in specific intervals of √ s from a lattice computation of the current correlator requires performing an inverse Laplace transform.Such lattice results will have uncertainties and be obtained at a finite number of points, making this inverse a notoriously illposed problem.
A number of interesting methods have been advocated for relating the two approaches.For instance, Ref. [22] proposes a modification of the Backus-Gilbert method that involves reconstructing, not the spectral function itself, but rather the spectral function convoluted, at individual values of √ s, with a Gaussian whose narrowness is limited by the statistical and systematic uncertainties on the lattice results.A first, application of this method was performed in Ref. [23].In Ref. [24] the authors propose a physics-and data-constrained dispersive representation of the pion electromagnetic form factor.Then, in Ref. [25] they use it to study the consequences of a value of a LO-HVP   µ   , such as the one obtained on the lattice in Ref. [6], on important observables impacted by this form factor. Another approach proposes to gain information about the R-ratio in different regions of √ s, via well-chosen linear combinations of so-called time-window observables computed on the lattice [26].Yet another approach proposes to use spectral-width sumrules to constrain R(s), in narrow regions of √ s, using the lattice correlator [27].
Here we propose to begin more modestly in terms of Rratio reconstruction.However, we do build a quantitative measure of comparison between the two approaches into our formalism and include it in our results.
To provide a relevant comparison between lattice and data-driven results, the lattice results used must have subpercent uncertainties.At such precisions, leadingorder strong-isospin-breaking and QED corrections are required.Moreover, if one is interested in probing the spectral function above the cc threshold, one must include the contributions of at least the four lightest quark flavors.For information above the b b threshold, a fifth quark flavour must be included.Of course, all continuum and infinite-volume limits must be taken in a controlled fashion.In addition, as the ranges of √ s that one is interested in reconstructing become narrower, the function with which the lattice results must be convoluted becomes more and more oscillatory.Thus, not only do the uncertainties have to be small, but also the statistical and systematic correlations between the lattice quantities have to be well known.
Here we use a limited form of the more general approach we propose, focusing on different intervals of √ s in the e + e − → hadrons spectrum, such as the ρ-peak region, the low-mass region below that peak, different high-mass regions, etc.The idea is to test the extent to which the lattice results are consistent with a modification of the spectrum that leads to a common rescaling of the observables of interest in the chosen region, and to determine the required rescaling factor.The simplest allowed modification would consist of directly rescaling the experimental R-ratio within the chosen region by that same factor.This would be a rather crude distortion of the experimental spectral function that can be viewed as a first approximation to a more physical modification.In fact, because the constraints that we consider pertain only to integrals of the R-ratio, there remains a significant amount of freedom in the shape of the corresponding modification.Therefore, many other modifications are possible, some of which may actually be physical.This point is discussed in more detail in Appendix E 1.
Beyond complete and precise lattice results for the HVP quantities discussed above, our proposed reconstruction and comparison approach also requires stateof-the-art determinations of the contributions to those quantities from different √ s-intervals in the data-driven approach, including a reliable quantification of correlations.For this we use the DHMZ methodology implemented in the HVPTools software [3,[70][71][72].
Another feature of our analysis is that we include "uncertainties on the uncertainties" for the lattice results.That is, we determine statistical and systematic uncertainties on the covariance matrix between the different lattice observables considered.This is important, because reconstruction methods are very sensitive to uncertainties and correlations in the input data, a corollary of the ill-posed nature of inverse methods.The same should eventually be done for the correlations between the corresponding data-driven observables.This is left for future work.
The remainder of the paper is organized as follows.Sec.II briefly presents the lattice correlation function and how it can be obtained from the R-ratio.The observables that we use in the comparison of the lattice and data-driven approaches are defined in Sec.III, including the window observables put forward in Ref. [8].In Sec.IV, we introduce the comparison methodology for testing the lattice correlation function with data-driven results.More importantly, we present our method for determining the size of the modification of the experimental R-ratio, in a given √ s-interval, that would be required to reconcile the lattice and data-driven values of the observables of interest.Sec.V presents results of applications of these methods to the comparison of the two approaches.This is followed by our conclusions in Sec.VI, in which we also summarize the main results of our study.
In addition, we provide a number of appendices.Appendix A presents the data-driven determination of the relevant observables and of their correlations and Appendix B, a summary of the lattice determination of these quantities, including uncertainties on the covariance matrix.In Appendix C we discuss, in more detail, the method used in this paper for determining the size of the possible modifications of the experimental Rratio and its possible pitfalls.Appendix D is dedicated to showing the stability of this method with respect to the averaging procedure which it calls upon.Appendix E presents various generalizations of this method, which can provide increasingly refined reconstructions of the R-ratio required to accommodate lattice results, as more of these results become available and/or their precision increases.Finally, in Appendix F, some of the limitations on the number of observables that can be studied in the data-driven approach are discussed.

II. THE HVP CORRELATOR FROM THE LATTICE AND THE R-RATIO
R(s) ≡ σ(e + e − (s) → hadrons(+γ)) 4πα 2 (s)/(3s) , where the denominator is the Born cross section for e + e − → µ + µ − in the massless limit.In that denominator, it is the electromagnetic coupling at c.o.m. energy squared, s, that is used, to remove unwanted vacuum polarization effects in the numerator.In particular, in comparisons of R(s) and C(t), this means that it is the one-photon-irreducible part of the latter that is relevant.We denote it C 1γI (t).Moreover, the cross section in the numerator of Eq. ( 2) includes the final-state radiation (FSR) of photons.While these photons imply that effects of higher order in α are included, they are kept to obtain an infrared-safe cross section at next-to-leading order (NLO) in α.Moreover, in applications, the contributions of these photons can consistently be taken into account in higher-order calculations.From R(s), it is straightforward to obtain C 1γI (t), after invoking the optical theorem that relates the e + e − → hadrons cross section to the spectral function associated with C 1γI (t).The result is proportional to the Laplace transform of √ sR(s) [21]: The integral diverges at short distances like t −3 , up to logarithms.However, in physical quantities, C 1γI (t) is multiplied by weights that vanish faster than t 3 as t → 0.

III. OBSERVABLES OF INTEREST
In principle, one could make a direct comparison of C 1γI (t) as a function of Euclidean time t, computed on the lattice and in the data-driven approach, via Eq.(3).However, C 1γI (t) has not yet been computed on the lattice, as a function of t, in the continuum and infinitevolume limits at the physical mass point.Here, to compare the two approaches, we focus on derived quantities, computed on the lattice with high precision, that are of direct phenomenological interest or that contribute to such quantities.
The LO-HVP contribution to the anomalous magnetic moment of the muon can be obtained from C 1γI (t) via [21]: where the kernel is (5) ω(r) = [r + 2 − r(r + 4)] 2 / r(r + 4) and α is the finestructure constant in the Thomson limit.In fact, it is a discretized version of Eq. ( 4) that is the basis of the lattice calculation of a LO-HVP µ .When the data-driven approach is used, a LO-HVP µ is computed directly from R(s) as an integral over s (see Appendix A).A comparison with lattice results gives a first means of confronting the two approaches.However, in the case of a disagreement, it is important to understand more precisely how different Euclidean times, t, in C 1γI (t) computed on the lattice, contribute to this tension.Indeed, different length scales on the lattice are subject to different statistical and systematic uncertainties, and knowing which ones agree or disagree with the data-driven approach may help point to aspects of the lattice calculations that require more attention.
To shed light on this issue, particularly useful quantities are the so-called Euclidean-time window contributions to a LO-HVP µ [8,26].These are restrictions, of the integral in Eq. ( 4), to time-intervals of interest, using a smoothed-out Heaviside function that limits edge effects when dealing with discrete times, as is the case on the lattice.This function is taken to be θ(t; ∆) where ∆ > 0 determines the timerange over which the function transitions from 0 to 1.The desired, restricted integrals are then defined as: where the window function, W (t; t i , t f , ∆), takes on different forms depending on the values of t i and t f : where we take ∆ = 0.15 fm, as advocaated in [8] Note that W (t; t i , t f , ∆) is defined such that a LO-HVP µ,win (0, ∞) = a LO-HVP µ .Also, in the following, we call a LO-HVP µ,win the Euclidean-time window contribution to a LO-HVP µ from the so-called "intermediate" time interval [0.4,1.0] fm, i.e.Eqs. ( 6) and (7) with t i = 0.4 fm and t f = 1.0 fm.Here we use a version of Eq. (5) in which the upper limit of the integral is Q 2 max = 3 GeV 2 for reasons explained in Ref. [19].
These time-window observables are particularly useful, because they are straightforward to compute in both the lattice and data-driven approaches.
On the lattice, their computation requires modifying the kernel in Eq. ( 4) used to compute a LO-HVP µ , that is using an appropriate discretization of Eq. (6).Of course, the difficulty of the continuum and infinite-volume extrapolations, as well as the quality of the statistical signal, will depend acutely on the chosen window.
With the R-ratio, one simply substitutes Eq. (3) into Eq.( 6), yielding: (9) A first quantitative comparison of the two approaches can proceed through the comparison of individual windows, as was done in Refs.[6,26].
In addition to time-window observables, one can further compare the lattice and data-driven approaches via the HVP function, Π(Q 2 ), at different values of the spacelike momentum q, with Q 2 = −q 2 ≥ 0 [19,20].Alternatively, one can consider the difference of this function at two values of Q 2 , to isolate an interval of interest.This has already be done in Ref. [6].Instead of directly considering Π(Q 2 ), which is convention dependent, here we focus on the hadronic running of the electromagnetic coupling in the on-shell scheme, to which it is related.Thus, we work here with the hadronic running of the fine structure constant in the five-flavor theory in the spacelike region, ∆ (5) had α(−Q 2 ).In terms of the lattice correlation function of Eq. ( 1), this quantity is given by and, in terms of the R-ratio, by So as to simplify notations in the following, we denote quantities computed on the lattice, a lat j , and their counterparts computed in the data-driven approach, a R j .Here the j indexes different moment integrals, such as the total a LO-HVP µ , window contributions to a LO-HVP µ , as well as the running of α between two spacelike virtualities.

IV. COMPARISON METHODOLOGY
In the following two subsections, we present in more detail the methodologies that we will use to compare and test the lattice and data-driven approaches in Sec.V.

A. Testing the lattice with R-ratio results
Using the notation introduced in the last paragraph of the previous section, the first step is to compare the observables of interest, a lat j and a R j , one by one, and to determine their level of compatibility.Depending on how localized these observables are in Euclidean time or in c.o.m. energy, one can isolate distance, virtuality or energy scales that may be problematic in each approach.
Since a more effective way of isolating possibly problematic regions of the experimental R-ratio is presented below in Sec.IV B, the observables a lat/R j that we focus on are chosen to emphasize features of the lattice correlator.Detailed results for the comparison of individual observables obtained from the lattice and the experimental R-ratio, are given in Sec.V A.
To go beyond a comparison of individual observables, we propose a more global comparison of the two approaches that includes not just one observable at a time, but many.Thus, we define the following χ 2 function: where j, k run over the observables considered in the comparison and where the a j are parameters.This function requires knowledge of the covariance matrices, C lat/R , between the observables obtained in each approach.Correlations between lattice and experimental R-ratio results are assumed to vanish, which is the case in the examples considered below.This more global measure of the compatibility of the two approaches is also important because it serves as a baseline for the improvements in the agreement brought by modifications of the experimental R-ratio considered in Sec.IV B.
For absolute uncertainties3 , the minimum of the χ 2 of Eq. ( 12), as a function of the parameters a i , is attained at and takes the minimum value, We take the p-value associated with the minimum value of χ 2 and the number of degrees of freedom (ndof) to be a measure of the overall agreement of the two approaches.Here, the ndof is simply the number of observables considered in the comparison.This measure of agreement makes sense only if C lat/R are well known, because the latter determine how much independent information is available in the lattice and R-ratio inputs and because the p-value obtained will depend strongly on them.This p-value may dilute some of the tensions observed in individual observables.However, in the presence of strong correlations, its value will better reflect the significance of having many observables disagree.Detailed results for this combined comparison of lattice and experimental R-ratio observables, via Eq.( 14), are given in Sec.V A.
Note that if the p-value corresponding to this minimum is acceptable, then the values of a j at the minimum of χ 2 correspond to a weighted average of the experimental R-ratio and lattice integrals.However, this averaging constrains the shape of the underlying spectral functions only in a very limited fashion.These aspects are discussed in more detail in Appendix E.

B. Testing the experimental R-ratio with lattice results: rescaling
While fully reconstructing the R-ratio from the HVP correlator computed on the lattice is an ill-posed problem, interesting information about R(s) can still be derived.
The method proposed here is applicable to any observable, a j , related to the HVP.It is meant to help isolate c.o.m. energy regions, in the experimental R-ratio, that may be responsible for tensions with the chosen lattice observables.Thus, we split up the observables obtained from the R-ratio, via e.g.Eqs. ( 8) and (11) Here the sum over b (i.e. over all of the intervals) covers the full support of R(s).In the case of a LO-HVP µ , of window observables and of the running of α between two spacelike momenta, the a R jb are obtained by restricting the integrals over s in Eqs. ( 8) and (11), respectively, to the interval I b .
When computing the Euclidean-time window observables, the edges of the corresponding intervals are smoothed out to account for the fact that, at finite lattice spacing, the density of points in C 1γI (t) is not that large.This is not necessary for R(s), because the density of measured points in s is significant, allowing us to consider sharp intervals I b , with a interpolation between the two points bracketing each boundary of the interval.
As discussed in the Introduction, in this paper we investigate whether a tension between lattice and datadriven results can be explained by a change in the experimental R-ratio that is consistent with a rescaling of the contributions, from one or more intervals in √ s, to the observables of interest.Thus, we wish to solve the following system of equations for γ b : taking into account all uncertainties and correlations among the lattice results, on the one hand, and the experimental R-ratio ones, on the other.If a given δ b ≡ (γ b −1) is set to zero a priori, the R-ratio integrals remain unchanged in the corresponding interval.
Here we consider the special case where a subset of the √ s-intervals, A, give contributions to the observables that are rescaled by a common factor γ while, for the complementary subset of intervals, B, the contributions are kept unchanged 4 .In such a case Eq. ( 16) becomes which implies, With more than one independent observable, a j , this becomes an over-constrained system, typically characterized by a χ 2 function: where C γ lat and C γ R are the covariance matrices of the γj rescaling coefficients, originating from the lattice and the dispersive uncertainties respectively.They are obtained through a linear propagation of uncertainties, from the lattice QCD and experimental R-ratio a j integrals, through Eq. (18).
The rescaling factor γ optimizing the constraints from Eq. ( 18) can be determined by minimizing the χ 2 function from Eq. ( 19) with respect to γ.This minimization yields an average of the input γj values, with weights given by the inverse of the sum of covariance matrices that appears in Eq. (19).Although very commonly used, this approach can yield biased results (see e.g.Ref. [73] and the discussion in Appendix C).In particular, such biases are caused by the uncertainties on the uncertainties and on their correlations, present both for the experimental measurements used in the dispersive approaches and for the lattice QCD results.In order to avoid these undesirable effects, we employ the weighted average of Eq. (C1), in which we set the off-diagonal elements of C γ lat and C γ R to zero, followed by a full propagation of the uncertainties with their correlations (see Appendix C for a detailed discussion).
We also consider a method based on a generalization of the χ 2 of Eq. ( 12), instead of the weighted averages discussed above.This method also allows us to study similar modifications of the experimental R-ratio.It is discussed in Appendix E 3.

C. Possible extensions of the methodology
In Appendix E, we consider comparisons and tests of the lattice and data-driven approaches that go beyond the rescaling of the experimental R-ratio integrals in a √ s-interval, by a single, common parameter γ.Moreover, in Appendix F we consider the possible advantages and limitations associated with the use of additional observables in the data-driven approach.

A. Testing the lattice with R-ratio results
As a first step, we compare individually a LO-HVP µ , the intermediate-time window observable a LO-HVP µ,win and δ(∆ (5) had α).The kernels for these quantities are shown as functions of t and √ s in Fig. 1.As is well known, the kernels for a LO-HVP µ is highly peaked for small values of √ s and large values of t.While the Euclidean-time kernel for a LO-HVP µ,win is strongly localized in t in the interval ranging approximately from 0.35 to 1.51 fm,5 its c.o.m. counterpart emphasizes a √ s region from threshold to around 3.2 GeV. 6Finally, the kernels for δ(∆ (5) had α) extends over all t, with an oscillation of period ∼ 1 fm, and over a large range of √ s, going from threshold to well beyond 5 GeV.The upshot is that the observables which we consider here probe the behavior of the lattice correlator C 1γI (t) at very different distance scales while the c.o.m. energies at which the R-ratio is probed have more overlap.Thus, the comparison of these observables tells us more about specific features of the lattice correlator than they do about those of the timelike spectrum.
The lattice [6] and data-driven results for the observables of interest are shown in Table I, together with their pairwise differences.While the preliminary lattice [6] and the data-driven results for δ(∆ (5) had α) are compatible within 1.4σ, there is a 2.0σ tension between the results for a LO-HVP µ .For a LO-HVP µ,win , the tension is much more significant, at 3.8σ.This number rises to 4.0σ when the lattice result of Ref. [6] is replaced by a weighted average of independent lattice determinations of this quantity [6,[9][10][11][12].All of these points are illustrated in Fig. 2.This large discrepancy points to a significant inconsistency between the approaches that must be resolved before giving a standard model prediction for a LO-HVP µ .Because the window observable probes intermediate Euclidean-time distances, its computation is particularly well suited to the lattice approach.The observable is significantly less noisy statistically than a LO-HVP µ , less sensitive to finite-volume effects, and less sensitive to shortdistance discretization effects, making its computation in lattice QCD particularly reliable.
In principle, the disagreement with the data-driven approach could be due to a common feature in the lattice calculations of Fig. 2 that would lead to an overlooked systematic error.For instance, one may worry that all the calculations are based on the time-momentum representation (TMR) of the current-current correlator given in Eq. (1).However, Ref. [75] explicitly checked that the connected contributions of the ud and s quarks to the window observable, obtained using a Lorentz-covariant coordinate-space representation (CCS), agree with those computed in the TMR to the level of 0.7%, in simulations with pion and kaon masses of ∼ 350 MeV and ∼ 450 MeV, respectively.Alternatively, one may question the universality of the continuum limit.However, as discussed in the caption of Fig. 2, the lattice results for [a LO-HVP µ,win ] ud iso , included in the average [6,[9][10][11][12] are obtained with different discretizations, thus testing this universality.Moreover, two of these analyses were blinded [11,12], strengthening this test.There are four additional calculations of the leading u and d quark contributions [8,[76][77][78].All nine calculations agree, as shown in the left panel of Fig. 2.
From the perspective of the R-ratio as a function of √ s, the range covered by the intermediate window is more spread out.Though it is not exclusively emphasized, this region includes the ρ peak, in which some disagreement between the different measurements is apparent [3,13].However, the difference between the KLOE [80][81][82] and the BABAR [83] measurements of the e + e − → π + π − (γ) cross section is taken into account as an additional systematic uncertainty, as described in Ref. [3].So, unless a large, unknown systematic uncertainty is present in the measurements of the e + e − → hadrons cross section [80][81][82][83], the strong tension with the lattice calculation of the intermediate window cannot be explained within the data-driven approach.The very recent measurement of the e + e − → π + π − (γ) cross section by the CMD-3 collaboration [13] helps bring the data-driven results for a LO-HVP µ and a LO-HVP µ,win in line with the lattice results for those quantities.However, it does so at the expense of a dramatic tension with previous measurements that is not yet understood.
As a preamble to the following discussion about the consequences, for the lattice correlator, of the tensions and agreement in the values of a LO-HVP µ , a LO-HVP µ,win and δ(∆ (5) had α) computed in the two approaches, it is important to remember that the current correlator is a smooth, monotonically decreasing function of t that behaves as t −3 for small t, up to logarithms, and falls off exponentially with exponent 2M π t at large t in infinite volume.
While the significance of the difference between the lattice and data-driven determinations of a LO-HVP µ is smaller than in a LO-HVP µ,win , the absolute difference itself is about twice as large.Thus, about half the discrepancy in a LO-HVP µ must come from the value of the lattice correlator for t below ∼ 0.35 fm and/or above ∼ 1.51 fm.An important caveat is that the uncertainty on this second half of the discrepancy is around 100% (see Table I), making this contribution to the discrepancy not significantly different from 0.
Combining this observation with the results for δ(∆ (5) had α), which agree in the two approaches, is more complicated.Indeed, the kernel of this observable is oscillatory in t and has support at all times.However, the kernel of δ(∆ (5) had α) has significant overlap with the one of a LO-HVP µ,win .This suggests that δ(∆ (5) had α) will receive an enhancement similar to the one of a LO-HVP µ,win , from the region delimited by the kernel of this window observable.Of course, that will not be true if all of the enhancement in a LO-HVP µ,win comes from times 1.2 fm < ∼ t < ∼ 1.5 fm, where the running-of-α kernel is suppressed.Assuming, for the moment, that the enhancement of δ(∆ (5) had α) is present, agreement on this quantity in the two approaches then requires that the lattice correlator be suppressed for some values of t outside the window.
This discussion leads to a situation where a LO-HVP µ suggests that the lattice correlator is enhanced for t below ∼ 0.35 fm and/or above ∼ 1.51 fm, while δ(∆ (5) had α) suggests the opposite.To try to reconcile the two statements it is useful to look at contributions to each of the two observables from the short-distance (SD) window (defined via the window function of Eq. ( 7), with t i = 0 fm and t f = 0.4 fm), the intermediate-distance (ID) one (denoted simply as window in this paper, and with t i = 0.4 fm and t f = 1 fm), and the long-distance (LD) one (with t = 1 fm to ∞).Using the data from Ref. [84] for a rough estimate of the ratio of the contributions of the SD:ID:LD windows, we find 10%:33%:57% for a LO-HVP µ and 70%:29%:1% for δ(∆ (5) had α).Thus, if the excess in δ(∆ (5) had α) from the ID window were to be compensated by a suppression of the correlator in the LD window alone, which contributes only 1%, this suppression would have to be significant.However, such a significant suppression would also reduce a LO-HVP µ , since the contribution of the LD window to that observable is dominant.To compensate that suppression would require a significant enhancement of the lattice correlator in the SD window, because that window only contributes at the level of 10% to a LO-HVP µ .But such an SD enhancement would make the lattice determination of δ(∆ (5) had α) significantly larger than that from the experimental R-ratio, given that the SD window is responsible for about 70% of the value of that observable.This discussion suggests that the SD window is not responsible for the additional enhancement in the lattice computation of a LO-HVP µ over that of a LO-HVP µ,win .It also suggests that an LD suppression of the lattice correlator cannot be responsible for mitigating a possible enhancement due to the ID window contribution to δ(∆ (5) had α).To summarize, we know for sure that the lattice correlator is enhanced, in the ID window region, compared to the one obtained from the data-driven approach.In addition, the arguments made above suggest that the lattice correlator is smaller in the SD window and larger in the LD one.
We now consider combined comparisons of the lattice and experimental R-ratio results for the above observables, via the minimization of the χ 2 function defined in Eq. (12).Beyond giving a more global measure of the compatibility of the two approaches, the resulting pvalues also provide baselines against which to compare those obtained when we study possible rescalings of integrals of the R-ratio in chosen √ s-intervals, in the following section.These comparisons require knowledge of the covariance matrices, among the observables considered, in each of the lattice and data-driven approaches.TABLE I. Comparison of lattice and data-driven results for the three main observables of interest in this paper.Each row contains the results and comparisons for the observable described in the first column.The second column provides the corresponding lattice result from Ref. [6] and the third, the data-driven one computed here, using the methods of Ref. [3].The fourth column displays the absolute difference between the two results, the fifth their relative difference, the sixth their difference in units of combined standard deviations and the last column, the corresponding p-value.The values in parentheses correspond to the total uncertainties of the corresponding quantities.Note that the powers of 10 in the observable column only apply to columns two through four.The three remaining columns have units specified in the corresponding column label.
Observable lattice [6]  a This result's continuum limit does not include the logarithmically-enhanced discretization uncertainties discovered subsequently in Ref. [79], nor was this quantity the focus of Ref. [6].However, we include it in the present study to illustrate how using a quantity which is complementary to a ] ud iso (green squares).For each group, only the most recent results are shown.BMW'20 [6], LM'20 [76], ABGP'22 [78] and FHM'23 [12] are obtained with different varieties of staggered fermions; Mainz'22 [9] with O(a)-improved Wilson, ETMC'22 [10] with twisted-mass and RBC/UKQCD'23 [11] with domain-wall fermions; χQCD'22 [77] with overlap valence quarks on either HISQ or domain-wall configurations.The filled squares correspond to fully independent results while the open ones are obtained using subsets of configurations from other calculations.LM'20 relies on a subset of the configurations used in FHM'23 and χQCD'22 on subsets of those used in FHM'23 and RBC/UKQCD'23.The green band corresponds to a weighted average of the fullyindependent (filled squares) lattice results for [a LO-HVP µ,win ] ud iso .The mean is performed without correlations in the determination of the weights while the uncertainty propagation assumes a 100% correlation between the total systematic error of the different calculations.The resulting correlated χ 2 /ndof is 2.3/4.Right panel: Comparison of the weighted average of lattice results (green filled circle and band) with a number of R-ratio results for the intermediate-window observable (red diamonds), including the one determined in this paper (filled red diamond and band).The weighted average for a LO-HVP µ,win is obtained from that for [a LO-HVP µ,win ] ud iso , to which we have added all other quark, QED and strong-isospin contributions from Ref. [6].The resulting average is a LO-HVP µ,win = (235.7 ± 0.8) × 10 −10 .The difference of this lattice average with the data-driven one of the present work is 4.0σ, shown as a horizontal arrow."WA" stands for "world average".
While the correlations between a LO-HVP µ and a LO-HVP µ,win , determined in lattice QCD, are obtained as described in Appendix B, the determination of the running of α is impacted by uncertainties that are, to good approximation, independent of those of the former.On the other hand, for the dispersive approach, the correlations among the uncertainties of the contributions to all moment integrals are derived as discussed in Appendix A.
We first consider a combined comparison of a LO-HVP µ and a LO-HVP µ,win .
The resulting χ 2 /ndof range from 14.4 +3.0 −2.1 /2 to 18.8 +2.0 −1.7 /2 when the two most extreme of the four systematic variations of the lattice covariance matrices, given in Appendix B 7 , are considered.The central values correspond to the nominal lattice sample.The uncertainties are statistical, obtained from the ±1σ quantile deviations from the median of the bootstrap distributions for the lattice covariance matrices 8 .Half the difference between the central values, i.e. 4.4/2 = 2.2, can be viewed as the lattice systematic uncertainty on these results.It is of comparable size to the statistical one.Unlike what has been just discussed for the lattice results, statistical and systematic uncertainties on the uncertainties associated with the data-driven approach are not yet included here, but should be in the future.
The corresponding p-values are 7 +10 −6 ×10 −4 and 8 +13 −6 × 10 −5 , respectively.These are very small and significantly smaller than the p-value for a LO-HVP µ alone in Table I.Nevertheless, the first of the above p-values can be as much as 17 times larger than the one in given in Table I for a LO-HVP µ,win and the second as much as a factor of 5 smaller, within one standard deviation.In any case, the probability that the lattice and data-driven results for a LO-HVP µ and a LO-HVP µ,win agree simultaneously is very small.
The disagreement reduces some when δ(∆ (5) had α) is added to the mix.In the lattice calculation, which is still preliminary, this quantity is negligibly correlated with a LO-HVP µ and a LO-HVP µ,win . While here the χ 2 remain essentially unchanged, ranging from 14.4 +3.0 −2.1 to 18.8 +2.0 −1.7 , the ndof increases from 2 to 3, leading to improved pvalues that range from 2.3 +4.0 −1.9 × 10 −3 and 3.0 had α) and, to a lesser extent, a LO-HVP µ is observed to dilute the disagreement between the lattice and data-driven approaches for a LO-HVP µ,win .This dilution is expected, because the disagreement between the two approaches is much smaller for both δ(∆ (5) had α) and a LO-HVP µ .However, in the presence of strong correlations, this combined measure will better reflect the significance of having agreement or disagreement among many observables.
Very similar results are obtained in the approach of Sec.IV B, when the γ rescaling coefficient in Eq. ( 16) is set to unity.This is not surprising because, with γ = 1, Eqs. ( 14) and (19) indicate that the χ 2 and ndof are the same in the two approaches, up to possible nonlinearities in the propagation of uncertainties.In particular, the minimal/maximal χ 2 /ndof values are stable, within 0.2 or less, with respect to the choice of the split in various √ s regions.When sampling the χ 2 /ndof values for 7 These two results correspond to the covariance matrices labelled 2 and 3, respectively. 8As discussed in Appendix B, it is important to note that the statistical covariance of the lattice data is determined from only 48 binned samples, which is smaller than the number of independent samples in the full dataset.However, as seen below, the statistical uncertainties on the lattice covariance matrix that arise due to this conservative binning are sufficiently small for the purposes of the present study.In future work, we will aim to improve this evaluation.
the various bootstrap replicas of the lattice covariance matrix, with either two or three moment integrals, the median of the resulting distribution is close to the nominal values (within 0.15 or less), while the variance and ±1 σ quantiles of the distribution are smaller than 1.5 and the ±2 σ quantiles are smaller than 3.0 (similar to the uncertainty values given earlier in this section).It is hence observed that the level of the tension between the dispersive-and lattice-based moment integrals is relatively stable with respect to the choice of the split in various √ s regions and/or statistical and systematic variations of the lattice covariance matrix.
Of course, other quantities related to hadronic vacuum polarization, that are computed in the lattice and datadriven approaches, can also be compared individually or added with correlations to the χ 2 , to further sharpen the comparison of the two approaches.

B. Testing the experimental R-ratio with lattice results: rescaling
In this section we present the results obtained when comparing the lattice QCD and dispersive results, employing the methodology discussed in Sec.IV B and Appendix C. With this methodology, we determine the size of modifications to the experimental R-ratio that lead to a rescaling of the contributions, from a given √ s-interval, to the observables of interest.As in the previous subsection, the study is performed for either two moment integrals (a LO-HVP µ and a LO-HVP µ,win ) or three of them (i.e.including, in addition, δ(∆ (5) had α)).Results are given for various scenarios, i.e.
√ s regions where the dispersive integrals are computed and the rescaling factors are determined.These include √ s regions below or above 0.63, 0.96, 1.1, 1.8 and 3.0 GeV respectively, as well as rescaling inside or outside the [0.63, 0.92] GeV interval.While the first set of intervals covers various low-and highmass regions, with different contributions to the dispersive integrals and different methodologies used for deriving them, the last one the region dominated by the ρ-resonance peak and the region complementary to it.The intervals are displayed in the left panel of Fig. 3.All uncertainties and correlations are accounted for, as discussed in the previous subsection.There we needed those associated with the moment integrals of the R-ratio over the full range of √ s.Here we need those corresponding to the contributions to moment integrals from the various √ s-intervals, whose determinations are also described in Appendix A.
Table II and Fig. 3 also present the results for the rescaling factor, γ 1 ≡ 1 + δ 1 , for the various √ s regions I 1 discussed in the previous paragraph, using either two or three moment integrals.The results are determined via the averages defined by the χ 2 function of Eq. ( 19), in which γ is replaced by γ 1 and with the weights of the input γj values obtained by setting all correlations to zero, for reasons discussed in Appendix D. Out of the four  II.The first panel displays, in dark grey, the various √ s-intervals considered here.The first row corresponds to the ρ-peak.In the second panel, we plot the p-values that indicate the compatibility of the rescaling hypothesis, in the given interval, with the lattice and data-driven results for a LO-HVP µ and a LO-HVP µ,win (2 observables, blue triangles) and with the results for the additional observable δ(∆ (5) had α) (3 observables, orange upside-down triangles).These p-values are obtained from the χ 2 defined in Eq. ( 19), using the appropriate ndof.The error bars are the statistical uncertainties resulting from those on the lattice covariance matrices determined in Appendix B (see also footnote 8).The corresponding uncertainty distributions are plotted in Figs. 4 and 6 (for the third and eleventh row of this plot respectively).The filled/empty points correspond to the largest/smallest p-value obtained from the four systematic variations of the lattice covariance matrix of Eqs.B31-B34.The vertical dashed line corresponds to 5%.We consider solutions with p-values above and close to that line to be compatible with the rescaling hypothesis in the corresponding interval, and those far below to be incompatible.The rightmost panel displays the rescaling percentages, δ1, corresponding to the p-values on the same row that have the same plot symbol.Here the uncertainties are the leading ones, given by the nominal values of the lattice and R-ratio covariance matrices, i.e. not including the small, additional uncertainty associated with the statistical uncertainty on the lattice covariance matrices.
systematic uncertainty variations of the lattice covariance matrices discussed in Appendix B 2, results are presented here for the two which induce the smallest/largest χ 2 values in the fit for each of the intervals.Indeed, depending on the fit configuration, the smallest χ 2 values are generally obtained for matrices "0" and "2", while the largest ones for matrix "3".
More specifically, Table II and Fig. 3 report the nominal results (obtained using the nominal lattice covariance matrices) for: the rescaling percentage δ 1 with its uncertainty propagated from the covariance matrices of the lattice QCD and dispersive results9 in Table II; the χ 2 /ndof obtained by injecting γ 1 into Eq.( 19); the corresponding p-value.
Using either two or three moment integrals, the rescaling of the integrals of the R-ratio by a common factor, in the √ s regions below 0.96, 1.1, 1.8 or 3.0 GeV, allows a good agreement between the lattice and R-ratio approaches.In those scenarios, the rescalings range from approximately 2 to 5%, the lower number corresponding to the larger intervals.In addition, a rescaling by around 5% in the ρ-peak interval, [0.63, 0.92] GeV, also restores the agreement.However, the tension persists when the rescaling is performed in the region below 0.63 GeV.
When using only the two observables, a LO-HVP µ and a LO-HVP µ,win , agreement is also restored through rescalings above either 0.63, 0.96, 1.1, 1.8 GeV, or outside the [0.63, 0.92] GeV interval, while tensions persist for a rescaling performed above 3.0 GeV.
Adding the third observable, δ(∆ had α), increases the constraints and tensions arise for a rescaling in any of the intervals that do not include the ρ-peak.
It is worth noting that the values of the rescaling percentage δ 1 for the various configurations presented in Table II are much larger than the (sub-)percent-level uncertainties of the experimental R-ratio in the corresponding √ s regions.Therefore, such shifts are not to be interpreted as a plausible way of solving the tension between lattice QCD and the dispersive approaches, but rather as a way of comparing various hypotheses for the possible source(s) of the tension.
Also given in Table II are the ±1σ statistical quantiles ) or three (adding δ(∆ had α)), as indicated in column 1.In these scenarios it is assumed that the integrals of the experimental R-ratio in the √ s-interval I1 are rescaled by a common γ1 ≡ 1 + δ1.This interval is given in column 2 and the nominal rescaling percentage δ1 (obtained through the weighted average discussed in Sec.IV B, using the nominal values of the lattice covariance matrices obtained in Appendix B), in column 4, with the corresponding nominal uncertainty propagated from the covariance matrices of the lattice QCD and dispersive results (C γ lat and C γ R ) indicated between "()".Column 3 indicates which of the four systematic uncertainty variations of the lattice covariance matrices give the best and worst fit qualities.Column 5 gives the nominal χ 2 /ndof for the scenario (obtained by injecting γ1 into Eq.( 19)) and column 6, the corresponding p-value.The values indicated between "[]" correspond to the ±1σ quantiles based on the lattice bootstrap replicas discussed in Appendix B (see also footnote 8).The last column indicates the shift on ∆α (5) had (M 2 Z ) induced by the δ1 rescaling.obtained for δ 1 , for the χ 2 /ndof and for the corresponding p-value, when considering their distributions obtained from the bootstrap replicas of the lattice covariance matrices computed in Appendix D (see also footnote 8).Examples of such distributions are shown in Figs.4-7.In the case of Fig. 3, which summarizes our results, those uncertainties are shown only for the p-values.

Number of observables
For δ 1 , in many scenarios the ±1σ quantiles based on the lattice bootstrap replicas are smaller than the primary uncertainty propagated from the covariance matrices of the lattice QCD and dispersive results, by about an order of magnitude or more.The changes associated with the four systematic variations of the lattice covariance matrix, are also small, though somewhat larger than the bootstrap ones.This is not necessarily surprising because they are uncertainties on the primary uncertainty.In fact, the latter are of only a few percent.Furthermore, the same statistical and systematic uncertainties on the χ 2 /ndof and on the corresponding p-values do not change the fact that a given fit can be considered "good" or "bad".These observations indicate that the statistical and systematic variations of the lattice covariance matrix do not change the conclusions qualitatively.
As mentioned in the previous subsection, the same bootstrap resampling can present some slight differences between the mean, median and nominal values of the corresponding distributions.In particular, the asymmetric, non-Gaussian tails (present especially for the p-values) impact the mean values (see e.g.Figs. 4 and 6).Still, these differences are well within the ±2σ quantiles (most often within the ±1σ quantiles) and do not impact the conclusions of the study.The asymmetric tails also induce some differences between the uncertainties obtained from asymmetric variances 10 , those from the ±1σ quantiles and those from the ±2σ ones (divided by 2).However, here also the conclusions of the study are not impacted by the slight differences among these three estimates of uncertainties.
In order to display the implications of observed normalization shifts for EW precision observable (EWPO) fits, in Table II we indicate the impact of the δ 1 rescaling on ∆α (5) had (M 2 Z ). 11As expected, the largest impacts on ∆α , when one also includes the 10 These variances are computed from the values in the distribution that are above and below the nominal one. 11While the uncertainty of the δ 1 rescaling factors is determined from the fits, we do not attempt a correlated propagation to (the uncertainty of) the corresponding shift of ∆α (5) had (M 2 Z ), nor a simple rescaling of the uncertainty of the latter.This follows the rescaling approach discussed in Ref. [18].There it was pointed out that, if such rescaling is necessary for the nominal values of the hadronic spectra, it is not obvious how the corresponding uncertainties should be computed.
hadronic contribution to the five-flavor running of the electromagnetic coupling, between −10 and −1 GeV 2 , the fit quality becomes very poor and the common rescaling of the contributions to the observables, in these high-√ s intervals, is strongly reduced too.Thus, in the presence of the δ(∆ (5) had α) constraint, the corresponding shifts on will not significantly impact the conclusions of EWPO fits.
As they are for the comparisons of Sec.V A, the conclusions of the studies performed here are stable against the statistical and systematic variations of the lattice covariance matrix.As shown in Appendix D, they are also stable with respect to the averaging methodology that is employed.Moreover, very similar results are obtained using the χ 2 approach to rescaling described in Appendix E 3.

VI. SUMMARY OF RESULTS AND CONCLUSIONS
We have presented a framework that enables a comparison of the primary ingredients that are used for determining HVP in the lattice QCD and data-driven approaches.The framework makes use of observables computable in both formalisms.These primary ingredients are: for the lattice approach, the zero three-momentum, current-current correlator of Eq. ( 1), studied as a function of Euclidean time, t; for the data-driven approach, the experimental R-ratio studied as a function of timelike c.o.m. energy, √ s (Eq.( 2)).The observables are integrals over t and √ s, weighted with appropriately chosen kernels.We alternatively call these observables or moment integrals.
The kernels can be chosen to be localized on the t-axis, as are the a LO-HVP µ time windows proposed in Ref. [8].Then by comparing the corresponding observables obtained in the two approaches, one can aim to isolate distance scales in the lattice correlator that may lie at the origin of the observed disagreement in the values, for instance, of a LO-HVP µ and a LO-HVP µ,win .Such comparisons are also interesting for HVP observables that have limited support in photon virtuality.These are typically related to the Adler function and the hadronic running of the electromagnetic coupling.They provide a complementary way to isolate relevant scales in the lattice correlator.More generally, our method applies to any observable related to HVP.
Our framework also allows us to gain information about the R-ratio in chosen regions of c.o.m. energy, from these same observables.The main challenge here comes from the fact that this requires solving a notoriously ill-posed, inverse problem, that arises because the lattice correlator is a Laplace transform of the R-ratio (Eq.( 2)).It also comes from the fact that observables which are localised in Euclidean time are much less so in c.o.m. energy.We solve the second problem by partitioning the √ s-axis from threshold to infinity into intervals,  4. Distributions of χ 2 (left) and p-value (right), obtained by sampling the bootstrap replicas of the lattice covariance matrix.These figures show the effects of different approaches to determining γ1 and its uncertainty in the low-mass region below 0.96 GeV.Two moment integrals are considered here, with the lattice covariance matrix "0".In the first row, we consider normalization fits using the full covariance matrices for determining the averaging weights-their comparison with the two plots in the following row, is discussed in Appendix D. In the middle row, we consider normalization fits with the averaging weights proportional to the inverse of the γj uncertainties squared, as considered in the present section.In the final row, we set the γ b normalization coefficients to unity, i.e. we perform no rescaling-these are discussed in Sec.V A. The blue horizontal error bar indicates the asymmetric variance of the distribution, computed on each side of the nominal value (blue point and blue dashed vertical line).The median (dashed-dotted black line), mean (black continuous line), 68.3% and 95.4% quantiles (coloured bands) of the distribution are also indicated.The number of degrees of freedom in the χ 2 calculation is indicated by the dotted pink line.Note that the scales in the last row of plots (corresponding to no rescaling) are very different from those of the two previous rows (in which rescaling is allowed) to account for the fact that the χ 2 and p-values with no rescaling are very poor.FIG. 5. Distributions of the normalization rescaling factor γ1 (left) and its uncertainty propagated from the covariance matrices of the lattice QCD and dispersive results (right), obtained by sampling the bootstrap replicas of the lattice covariance matrix.These figures show the effects of different approaches to determining γ1 and its uncertainty in the low-mass region below 0.96 GeV.Two moment integrals are considered here, with the lattice covariance matrix "0".In the first row, we consider normalization fits using the full covariance matrices of the γj for determining the averaging weights.In the second row, we consider normalization fits with the averaging weights proportional to the inverse of the γj uncertainties squared.The blue horizontal line indicates the asymmetric variance of the distribution, computed on each side of the nominal value (dashed blue vertical line).The median (dashed-dotted black line), mean (black continuous line), 68.3% and 95.4% quantiles (coloured bands) of the distribution are also indicated.and computing the contributions of each such interval to the observables of interest.The first challenge -the reconstruction of a continuous function from a finite number of its weighted integrals does not have, a priori, a unique solution -is addressed by limiting the information that we try to extract while monitoring the possible pitfalls associated with such a reconstruction.
In the present work, we focus on simple scenarios.We consider three observables related to HVP: a LO-HVP µ , a LO-HVP µ,win and the hadronic running of the electromagnetic coupling between spacelike virtualities 1 GeV 2 and 10 GeV 2 , δ(∆ (5) had α).This choice is motivated by the fact that these quantities have been determined on the lattice to the sub-percent level, with all relevant corrections [6].It is important to note that the result for δ(∆ (5) had α) is not the main focus of Ref. [6] and should be understood as preliminary.Nevertheless, we include it because it im-poses complementary constraints on the spectral function and demonstrates the flexibility of our approach.
Then, we use these observables to: I. isolate regions of t in which the lattice correlator leads to values of the observables that are inconsistent with those obtained in the data-driven approach; II. determine intervals of √ s in which a change in the experimental R-ratio could explain the observed discrepancies between the lattice and R-ratio approaches, by allowing a common rescaling of the the contributions from these intervals to the observables of interest.
Because our approach is based on the minimization of χ 2 functions, it allows us to quantify the compatibility of the observables in the two approaches, as well as FIG. 6. Distributions of χ 2 (left) and p-value (right), obtained by sampling the bootstrap replicas of the lattice covariance matrix.These figures show the effects of different approaches to determining γ1 and its uncertainty in the high-mass region above 3.0 GeV.Two moment integrals are considered here, with the lattice covariance matrix "3".In the first row, we consider normalization fits using the full covariance matrices for determining the averaging weights-their comparison with the two plots in the following row, is discussed in Appendix D. In the middle row, we consider normalization fits with the averaging weights proportional to the inverse of the γj uncertainties squared, as considered in the present section.In the final row, we set the γ b normalization coefficients to unity, i.e. we perform no rescaling-these are discussed in Sec.V A. The blue horizontal error bar indicates the asymmetric variance of the distribution, computed on each side of the nominal value (blue point and blue dashed vertical line).The median (dashed-dotted black line), mean (black continuous line), 68.3% and 95.4% quantiles (coloured bands) of the distribution are also indicated.The number of degrees of freedom in the χ 2 calculation is indicated by the dotted pink line.Note that the scales in the last row of plots (corresponding to no rescaling) are very different from those of the two previous rows (in which rescaling is allowed) to account for the fact that the χ 2 and p-values with no rescaling are very poor.7. Distributions of the normalization rescaling factor γ1 (left) and its uncertainty propagated from the covariance matrices of the lattice QCD and dispersive results (right), obtained by sampling the bootstrap replicas of the lattice covariance matrix.These figures show the effects of different approaches to determining γ1 and its uncertainty in the high-mass region above 3.0 GeV.Two moment integrals are considered here, with the lattice covariance matrix "3".In the first row, we consider normalization fits using the full covariance matrices of the γj for determining the averaging weights.In the second row, we consider normalization fits with the averaging weights proportional to the inverse of the γj uncertainties squared.The blue horizontal line indicates the asymmetric variance of the distribution, computed on each side of the nominal value (dashed blue vertical line).The median (dashed-dotted black line), mean (black continuous line), 68.3% and 95.4% quantiles (coloured bands) of the distribution are also indicated.
the consistency of the observables with our hypothesis that a change of the experimental R-ratio, in a particular √ s-interval, could explain the observed discrepancies.Clearly, the very simple model of item II tells us little about the shape of the deformations of the timelike spectral function that may be required.Moreover, it does not incorporate, a priori, physical constraints, such as those imposed by chiral perturbation theory at small c.o.m. energies, by perturbative QCD at large ones, or by analyticity and positivity at all energies.However, it is not necessarily incompatible with some of them either.In addition, it can help isolate √ s-intervals in which the experimental R-ratio may have to be modified to explain the lattice results.
Within a given approach, the observables obtained are not independent.Thus, it is critical to take into account correlations between them.This is necessary to correctly quantify the compatibility of the lattice and data-driven observables with one another.It is also particularly important when considering the inverse problem, because its solutions are very sensitive to these correlations.In fact, because of the exponential enhancement of fluctuations, it is also important to consider the stability of the results against uncertainties on the corresponding covariance matrices.We take into account all uncertainties and correlations and, for the lattice results, we also consider the effect of uncertainties on the corresponding covariance matrices.The consideration of such an effect, in the data-driven approach, is left for future work.
Using the methods described in Sec.IV A we find the results presented and discussed in Sec.V A. They represent our attempt to isolate regions in the lattice correlator, on the Euclidean time axis, that do not agree with the data-driven approach.
As already observed in Ref. [6], the lattice result for a LO-HVP µ is 2.1σ above the one obtained from the data-driven approach (corresponding to a p-value of 5 × 10 −2 ), the one for a LO-HVP µ,win , 3.8σ above (p-value of 1 × 10 −4 ), while the preliminary lattice result for δ(∆ (5) had α) is only 1.4σ above (p-value of 0.15).Moreover, the lattice result for the leading, isospin-limit, up-and-down quark contribution to the window observable has been independently confirmed by nine collaborations [6,[9][10][11][12], of which two also confirm the results for other contributions [9,11].This increases the tensions with the data-driven results to a conservative 4σ.
The discrepancies discussed in the previous paragraph are confirmed when correlations are taken into account via Eq.(12).Considering first a LO-HVP µ and a LO-HVP µ,win together we find p-values which range approximately from 2 × 10 −5 to 1.7 × 10 −3 , taking into account the maximum systematic variation and 1σ statistical fluctuation in the lattice covariance matrices.These are substantially smaller than the one for a LO-HVP µ alone.Adding to the comparison the preliminary lattice result for δ(∆ (5) had α), and its data-driven counterpart, reduces the overall tension to some degree.The p-values now range from approximately 10 −4 to 6 × 10 −3 , still making an overall agreement very unlikely.
Combing all of this information with the shape of the observable's kernels, in Sec.V A we argue that our results point to an excess of the lattice correlator in the intermediate-time range from 0.4 to 1.5 fm, a probable excess at time separations above 1.5 fm and a possible suppression in the short-time range below 0.4 fm.
Then, we turn to the implications that the comparison of the two approaches may have on the experimental Rratio.Again, the tests considered take into account all systematic errors and correlations, as well as uncertainties on uncertainties in the lattice results.
Using the averaging procedure of Sec.IV B, we find that the lattice predictions, for a LO-HVP µ and the corresponding a LO-HVP µ,win (but excluding δ(∆ (5) had α)), can be reproduced within uncertainties via a modification of the experimental R-ratio in various intervals of √ s.More specifically, if the lower end of these intervals begins at threshold, the intervals must include the ρ peak.Indeed, for intervals whose upper bound is anywhere between √ s = 0.96 to 3.0 GeV, the modifications have to be such that they increase the contributions to both the data-driven a LO-HVP µ and a LO-HVP µ,win in those intervals by approximately 2 to 5%, the lower value corresponding to the larger intervals.These conclusions allow scenarios, such as those advocated in [85], in which the R-ratio would have to receive modifications above 1 GeV to explain the lattice results for a LO-HVP µ and a LO-HVP µ,win .Interestingly, a modification of only the ρ-peak region is also compatible with the lattice results, but requires around a 5% increase of the contributions to the two observables from that region.It is worth noting that such an increase is comparable to the one claimed by CMD-3 in a very similar region [13].However, without having worked out the correlations between CMD-3 and the other e + e − → hadrons measurements, it is impossible to make a more quantitative statement.Moreover, in the presence of large, unexplained discrepancies between measurements, and thus of large unknown systematics, one cannot claim to control such correlations and uncertainties.
Somewhat surprisingly, a scenario in which the experimental R-ratio is modified in an interval from threshold to infinity from which the ρ peak is excluded is also permitted, with p-values above ∼ 10% in the two-observable case-i.e.beginning right below the cc resonance region.Moreover, modifications of the experimental R-ratio, in intervals that extend to infinity and begin at √ s as small as 0.63 GeV and as large as 1.8 GeV, are also allowed.In that case, the contributions to the data-driven a LO-HVP µ and a LO-HVP µ,win must be increased anywhere between approximately 2% to 40%.On the other hand, an explanation of the discrepancy between the two approaches via a modification of the experimental R-ratio only above 3 GeV is excluded with high probability.
In order to see if some of the above scenarios can be eliminated, we consider adding the observable δ(∆ (5) had α).This observable has a footprint along the √ s-axis that is significantly more extended towards larger c.o.m. energies than those of a LO-HVP µ and a LO-HVP µ,win .We find that the addition of this observable eliminates, with significant probability, all scenarios in which the modification of the experimental R-ratio excludes the ρ peak.Nevertheless, the scenario including only the ρ-peak interval, and the one in which the interval begins below the ρ peak and extends to infinity, both survive this additional constraint.Moreover, these possible modifications of the experimental R-ratio imply a rescaling of the data-driven results, for all three observables, by a common factor between about 1.02 and 1.06.In particular, the scenario including only the ρ peak is barely changed by the addition of the constraint from δ(∆ (5) had α), making it a robust explanation for the discrepancy between the lattice and data-driven approaches.
Not surprisingly, very similar conclusions are obtained using the fitted-moment formalism given in Appendix E 3. There, however, the natural interpretation of the results is slightly different.In that approach, the contributions of different √ s-intervals to the various observables become fit parameters, as are the rescaling factors.Thus, under the assumption that the proposed rescaling scenario is correct, as should also be verified by the corresponding χ 2 /ndof, the fitted contributions to the observables can be understood as a combination of lattice and data-driven results for these contributions.
It is important to remember that all of the modifications to the experimental R-ratio discussed above are not consistent with the uncertainties of the timelike spectrum.This is made obvious by the very small p-values, discussed above, that correspond to the minimum of the χ 2 of Eq. ( 12), which measures the combined compatibility, of the lattice and data-driven results, for a LO-HVP µ , a LO-HVP µ,win and δ(∆ (5) had α).As discussed in Appendix E, the methods presented above can be generalized to include more and more information about HVP from the lattice and data-driven approaches.This additional information will allow us to refine our understanding of the discrepancies between the two approaches, both in terms of the Euclidean correlation function and the experimental R-ratio.Our methods can also be generalized to include more theoretical constraints, such as those discussed earlier in this section.Moreover, the methods presented in Appendix E supply means of obtaining combinations of lattice and data-driven results for various quantities of interest related to HVP.These combinations will only be as valid as the model used to reconcile the two approaches is well motivated.
Finally, in Appendix F we argue that the number of linearly independent moments that can be obtained with reasonable precision, from the experimental R-ratio, is less than 10 for the types of observables considered here, which should be computable with sub-percent precision in the lattice approach.
Appendix A: Data-driven determination of the observables and their correlations Computing the contributions to the moments integrals from the low energy region, below 1.8 GeV, requires combining the spectra measured by various experiments.These measurements are performed either through a scan of the energy in the centre-of-mass of the collider, as is the case of e.g. the CMD-2 and SND experiments, or using the ISR technique, in the case of KLOE and BABAR (see e.g.Refs.[14,[80][81][82][83][86][87][88] for the π + π − channel, dominant in the case of a LO-HVP µ ).In the latter case, selecting events with a hard photon emitted from the initial state allows such studies to cover a broad range of invariant masses for the hadronic final state system.Actually, BABAR covers the full energy range of interest, due to selecting events with the hard photon reconstructed in the detector and to the high energy in the centre-of-mass of the collider.While the KLOE ′ 08 [80] and KLOE ′ 12 [82] measurements consider events with the hard ISR photon along the beam, the KLOE ′ 10 measurement uses large-angle ISR photons [81].Furthermore, in measurements where the ISR luminosity is evaluated in-situ based on e + e − → µ + µ − events [82,83,88], several systematic uncertainties are reduced.Still, among the currently available results, only the BABAR study uses the reconstructed extra photons (performing a measurement reconstructed at "NLO" and inclusive of any number of additional photons), necessary at the sub-percent level precision aimed for here.The differences between the sampled event topologies, the detector technologies, and especially the analysis methodologies, result in systematic tensions between the BABAR and KLOE measurements [3,72], but also between the KLOE measurements themselves (see Ref. [5] and references therein).Since these tensions are still to be understood and fixed, they require special care in the current combinations.
The channel-by-channel combinations in the DHMZ approach (implemented in the HVPTools software) [3,[70][71][72] are performed in fine energy ranges, using a splinebased procedure.It takes into account the information on the uncertainties and their correlations between bins/points and between experiments.The DHMZ procedure also accounts for correlations between the various different exclusive channels, when summing their contributions, which induces a necessary enhancement of the total uncertainty.Furthermore, the local tensions between the measurements are taken into account based on the computation of a χ 2 /ndof in fine energy ranges, used for the rescaling of the uncertainties.However, in the presence of systematic tensions, as the ones between the BABAR and the KLOE measurements in the ππ channel, the local rescaling of the uncertainties is not sufficient.For this reason, in the most recent DHMZ study [3], these tensions are taken into account by treating half of the difference between integrals without either BABAR or KLOE as an extra uncertainty.This yields the dominant uncertainty on a LO-HVP µ , amounting to 2.8 • 10 −10 .The DHMZ approaches for accounting for correlations between channels and for the BABAR-KLOE tension in the ππ channel, have also been adopted for the result of the Theory Initiative White Paper [5].
The procedure described above is used to derive the experimental values of R e + e − →hadrons (s), together with its covariance matrix, up to 1.8 GeV.It is completed with contributions from inclusive experimental measurements in the range between 3.7 and 5 GeV, combined using the same methodology as above, with contributions from narrow resonances, as well as from perturbative QCD for the continuum [3].An uncertainty accounting for possible duality violation effects is also included.A direct integration of the resulting R e + e − →hadrons (s) on the full mass range, together with a linear propagation of the uncertainties with their correlations, yields the moment integrals with the kernels of interest here, together with their covariance matrix.

Appendix B: Determination of lattice observables and of their correlations
The lattice QCD results used in this work were obtained in Ref. [6].In order to calculate the two-point function C(t) defined in Eq. ( 1), we need to evaluate the current correlator ⟨J µ (x)J μ(x)⟩ in our lattice regularization.We define the regularized current correlator as the second differential of the partition function with respect to an external photon field Thus, we use the conserved current at source and sink.The correlator can be evaluated as plus a contact term that gives no contribution to the observables of interest here.C conn µ,x,μ,x and C disc µ,x,μ,x arise respectively from the connected and disconnected Wick contractions of the quark fields and are given by where Tr is the spin-color trace, P x is a projection operator that sets all sites other than x to zero, M −1 f is the propagator for a quark of mass m f on a gauge background V U e ieq f A , and These calculations are performed gauge-ensemble by gauge-ensemble.The SU(3) configurations are generated at the isospin-symmetric point, where m u = m d and e = 0. Leading-order QED and strong-isospin-breaking (SIB) corrections are then added as perturbations.
The leading QED corrections are obtained as follows.For each SU(3) gauge configuration, a stochastic U (1) photon field is generated with the non-compact, QED L action.The U(1) fields are exponentiated into link variables that we multiply with the corresponding SU(3) links.The covariant Dirac operator is then inverted on each of the corresponding gauge configurations and the resulting quark propagators are appropriately contracted.The first and second derivatives of the contractions, with respect to the unit of electric charge e, are obtained by computing the contractions for +e and −e and by taking the appropriate finite differences.The expressions for the first and second derivatives of the fermion determinant, in terms of quark propagators, are computed analytically.For the second derivative, around 2000 photon fields are used.
These first and second derivatives are appropriately combined, with the e = 0 contractions as well, to obtain all O(e 2 ) contractions relevant for the current correlator.The desired correlators are then obtained by averaging the various contractions over the SU(3) and U(1) gauge configurations of each ensemble.
At leading order, the fermion determinant does not contribute SIB corrections.The valence contributions are obtained by computing the first derivative of the contractions, with respect to light-quark mass, and then multiplying it by the u-and-d quark mass difference.As above, the desired correlators are then obtained by averaging the various contractions over the SU(3) and U(1) gauge configurations of each ensemble.
In order to obtain physically relevant values, we perform global fits to the lattice spacing, quark mass, and electric charge dependence of a LO-HVP µ and a LO-HVP µ,win j.These fits take the form where the first line corresponds to a definition of the isospin-symmetric limit, and the second line gives the isospin breaking effects.The X l and X s variables describe the deviation from the physical light and strange mass where "ph" indicates the physical value.The remaining X variables measure the distance to the isospinsymmetric limit where e v and e s are the valence and sea electric charges, respectively.The meson masses are defined as The coefficients A, B, . . . in Eq. (B6) are specific to the observable Y .They can depend on the lattice spacing, and also on the X variables defined above, in particular we use: For any ensemble, a is given by the ratio of the Ω mass, in lattice units, to its experimental value.For the A coefficient of observables which give a large contribution to the final result, we include powers n of the strong coupling α s (1/a) in the lattice-spacing dependence [89].We take both the commonly used n = 0 (corresponding to the usual polynomial expansion in a 2 ), and n = 3, as suggested in Ref. [90].For the strong coupling, we use its four-flavor, MS value at scale 1/a.We determine this value from the world average value of α s (M Z ) [91], by running down from M Z to 1/a in five-loop perturbation theory [92], taking into account four-loop threshold corrections [93] at the b-quark mass that are given in Ref. [91].
The choice of n is included as a source of systematic uncertainty in the analysis described below.In addition, we vary the polynomial order of the coefficients as an additional source of uncertainty.As required by the data, we take the A coefficient to be linear, quadratic, or, in some cases cubic in a 2 α s (1/a) n .For all other coefficients, the lattice spacing dependence is assumed to be a function of a 2 only.In the D and E coefficients, up to quadratic dependencies are used, in all other cases only a linear one is needed.Depending on the fit qualities, some of these parameters will be set to zero.
We perform simulations at six different values of the lattice spacing, at a range of quark masses scattered around the isospin-symmetric point.We exclude zero or more of the coarsest lattice spacings from the analysis, and the resulting variation in the global fit provides another source of systematic uncertainty.

Estimating the covariance matrix
The results of this lattice calculation are subject to both statistical uncertainties arising from the Monte Carlo sampling of the QCD gauge fields, and systematic uncertainties arising from choices made during the analysis.We estimate the variance and covariance arising from both types of uncertainty by using a combination of statistical resampling, and histogramming of systematic variations [6,94,95].
The statistical covariance matrix is computed through jackknife resampling of the two-point, electromagnetic current-current correlator C(t).The gauge field configurations are computed through a Markov chain algorithm that induces an auto-correlation between subsequent configurations.We suppress this effect by introducing a blocking procedure, where the individual configurations are grouped together into N J blocks and then the average of C(t) on each block is used in the resampling procedure instead of the value on each configuration.To simplify the analysis, we take N J to be the same on all ensembles.Specifically, we take N J = 48 which gives typical blocks of 100 or more configurations, larger than the autocorrelation time of the topological charge (which is ∼ 20 on the finest ensembles).
Systematic uncertainties arise from a number of choices that must be made throughout the analysis.For example, the choice of appropriate fit ranges, the experimental value of the Ω-baryon mass used to set the scale, or the particular fit form used to perform the global fit.We call the global fit resulting from a particular set of these choices an analysis.In order to quantify the systematic uncertainties, we need to estimate the variance of the final value across all analyses arising from reasonable choices.To do this we follow the strategy used in Ref. [6].
As the full value of a LO-HVP µ and the window observable a LO-HVP µ,win are evaluated with similar techniques on the same set of gauge configurations, a number of choices in the analyses for a LO-HVP µ and a LO-HVP µ,win are similar and hence introduce correlations.Therefore, the procedure described in Ref. [6] had to be adapted to take into account these correlations.δ(∆ (5) had α) is assumed to have negligible correlations with a LO-HVP µ and a LO-HVP µ,win because it is mostly determined by contributions from shorter distances.
The systematic uncertainties affecting any given analysis can be arranged into two groups: those that induce correlations, and those that induce no correlation.The systematics that induce no correlation are further subdivided into those that enter with a flat weight, and those that are weighted with the Akaike Information Criterion: where the χ 2 , the number of fit parameters n par and the number of data points n data describe the global fit corresponding to the given analysis.These relative AIC weights are employed between analyses differing in the choice of fit functions or in the lattice spacings included.For other systematics, including all systematics that induce correlations, a flat weight is employed.The choices that are assumed to introduce correlations are: the fit ranges in the extraction of masses from correlation functions, the experimental input for the Ω mass, the taste improvement time-ranges, and the choice of t c in the bounding method for the extraction of a LO-HVP µ .Details of this procedure and explanation of the individual contributions to the systematic uncertainties can be found in Ref. [6].In the following, we denote with A(ϕ, ψ flat , ψ AIC ) the result for some observable A, either a LO-HVP µ or a LO-HVP µ,win , with a specific set of systematics.Here, ϕ is a label for a set of specific choices that will induce correlations between a LO-HVP µ and a LO-HVP µ,win and the labels ψ flat and ψ AIC refer to choices that are made independently and are weighted either by a flat weight or by an AIC weight.In this notation, the histogram for a single observable A can be expressed as where N 1 (x, µ, σ) is the probability density function of the normal distribution with mean µ and standard deviation σ.A(ϕ, ψ flat , ψ AIC ) refers to the mean value of the observable A in the analysis with the systematic variation indicated in parentheses.In the same sense, σ A(ϕ,ψ flat ,ψ AIC ) corresponds to the statistical uncertainty on this analysis determined from the jackknife samples.The weight function w(ϕ, ψ flat , ψ AIC ) = AIC(ϕ, ψ flat , ψ AIC ) ensures a correct relative weight between the analyses.
Here, AIC(• • • ) denotes the AIC weight of the indicated analysis.This procedure can be generalized to the study of the correlations between a LO-HVP µ and a LO-HVP µ,win by introducing a two-dimensional histogram: Here, C stat refers to the statistical covariance matrix determined from the jackknife samples for each individual combination of analyses and N 2 is the bivariate normal distribution with this covariance matrix.A and B refer , or individual flavor and contraction contributions to them.
We treated the systematic error introduced by the taste correction as a special case-both a LO-HVP µ and a LO-HVP µ,win are affected by this systematic in a correlated way.For a LO-HVP µ , four time-ranges for taste improvements to the correlation function are used, starting at the Euclidean times 0.4 fm, 0.7 fm, 1.0 fm, and 1.3 fm.But since the window observable a LO-HVP µ,win is not affected by the taste improvement in the last range, for a LO-HVP µ,win the two last ranges give the same results.Therefore, using all four ranges for both observables would bias the window observable.To avoid this but still retain a measure of correlation, we considered 12 systematic combinations where, for each matrix, we fixed the range of the taste correction to a LO-HVP µ and a LO-HVP µ,win .The starting times for the taste correction ranges in all 12 cases are shown in Table III.Then, we averaged the 12 covariance matrices.
As there is a large number of choices with many of these inducing a correlation, explicitly evaluating Eq. (B25), as a function of a LO-HVP .From these three histograms, we determined the median and the interquantile distances d σ X , d σ Y and d σ Z .Here, σ indicates the number of standard deviations the interquantile distances would correspond to, in the case of normal distributions.Using these results, we construct the covariance matrices as To avoid any bias from the details of the histogram bin-ning, we considered a large number of 2000 bins per histogram.The upper and lower bounds of the histogram were chosen such that they do not influence the results.These values can be found in Table IV.The resulting histograms for X and Y are plotted as black curves on the upper and on the right sides of the heat maps of Fig. 8.For illustrative purposes, we also calculated twodimensional histograms with a smaller number of 1000 bins in each direction.These two-dimensional histograms are also shown Fig. 8.The bins were filled by sampling 5 random points for each summand in Eq. (B25).Note that the random sampling is only included in the visualization and has no impact on the end result.It is noteworthy that in the case of the light quark contribution to the window quantity, a double peak structure appears.This structure is due to the variation of the parameter n = 0, 3 in Eq. (B15).The effect of the non-gaussianity of the distribution introduced by this structure is taken into account by looking at both the σ and 2σ quantiles.
In order to separate the systematic uncertainties and correlations from the statistical ones, we repeat, in analogy to the procedure in Ref. [6], the full procedure while replacing C with λC in Eq. (B25).We choose λ = 2.Then, we solve the equations , extrapolated to a = 0 and interpolated to the physical mass point.

Systematic and statistical uncertainties on the covariance matrix
The specific value taken by the estimated covariance matrix can have an important effect on the results described in the main paper, in particular on those for χ 2 , for the corresponding p-values, as well as those resulting from the inverse problem.As such, we found it important to quantify the uncertainties in this estimate.We identified two major sources of systematic uncertainty in our covariance estimate: non-normal distributions, and correlation of continuum extrapolations.Each of these sources is described below, along with our approach to estimating the associated uncertainty.
Our jackknife estimate of the statistical covariance matrix for each analysis has itself a statistical uncertainty.We estimate this uncertainty by applying an approximate second-order resampling technique.We observe that the jackknife bias in a LO-HVP µ and a LO-HVP µ,win is small.Hence, we can approximate the true bootstrap distribution by bootstrapping the reconstructed samples X i = N J X − (N J − 1) Xi , where X is the ensemble mean and Xi is the ith jackknife.We construct N B = 1000 bootstraps in this way, and then on each bootstrap we .These are obtained via the specialization of Eq. (B25) to these individual contributions, as described in the text.Each of the three plots has a heat map in its center, illustrating the two-dimensional covariance distribution of the total, combined uncertainty.As usual, colors from white to blue represent higher to lower probability densities.These densities are calculated using the random sampling discussed in the text.The solid and dashed ellipses in these maps indicate the one and two sigma errors determined from the one and two sigma quantiles.In the upper and right panels around the heat maps, histograms of the marginalized distributions are shown.These are calculated exactly without random sampling.The two green bands show the one and two sigma uncertainties determined by the one and two sigma quantiles.compute the jackknife statistical covariance and perform the histogramming technique described above to obtain the full covariance matrix.In this way we compute the bootstrapped distribution of our final covariance matrix.
To obtain a covariance from the histogram, we make the assumption that the histogram is normally distributed.Under this assumption, the one-sigma quantiles described above will fall one standard deviation below and above the median.Any deviation away from the normal distribution will hence introduce a systematic uncertainty.One example of such a non-normal distribution is the double-peak structure observed in Fig. 8a.This double peak arises from the two-point systematic related to the power n of the strong coupling in the description of the lattice spacing dependence of a LO-HVP µ,win,l (i.e. the connected light (ud) quark contribution to the window observable) given in Eq. (B15).It is an artefact of how we sample the errors, and not intrinsic to the underlying distribution, but it still introduces an uncertainty in how the covariance is defined.To quantify this uncertainty, we also consider the covariance matrix obtained by considering the two-sigma quantiles to be two standard deviations below and above the median.If the distribution were normal, the covariance matrix obtained in each of these ways would be identical.We take the difference between these two determinations of the covariance matrix as a two-point systematic.
In producing the histograms we allow the form of the continuum fit to vary independently for a LO-HVP µ and a LO-HVP µ,win , as there is no a priori reason that the fit forms should be the same.Hence, by construction we obtain no correlation between systematic uncertainties arising from the continuum limit of each quantity.However, since a LO-HVP µ can be considered to be the sum of W = a LO-HVP µ,win and its complement C, it is reasonable to surmise that there may be some correlation between their continuum limits.To estimate this potential correlation, we consider taking the continuum limits of W and C independently, and adding them together to obtain T = W + C = a µ .In this scenario, let W have variance (dW ) 2 c associated with the continuum limit, and similarly for C. Then the variance of T is where (dT ) 2 o is the variance of T due to uncertainty sources other than the continuum extrapolation.Under these assumptions, the covariance matrix between T and W is where C o T W describes the covariance not associated with the continuum extrapolation.In practice, since the continuum limit of the combined T = W + C is taken instead, the C component may interfere and partially spoil the correlation.As such we expect the true covariance to lie somewhere between this case, and the case where this contribution is neglected.We take as an additional twopoint systematic covariance matrices including either 0% or 100% of this estimated continuum limit contribution.
The connected contribution from valence charm quarks was computed in Ref. [96], and the data was not available to perform the covariance analysis described here.To address this, we assume that the correlation between a LO-HVP µ,c and a LO-HVP µ,c,win lies somewhere between 0% and 100%, and treat the difference between these two extremes as an additional uncertainty.The charm contribution is a very small part of the total variance so this uncertainty is small.It is closely related to the uncertainty from the correlation of continuum extrapolations and as such is treated as part of the same two-point systematic.
Finally, the finite-size effects computed in Ref. [6] have a large contribution from the continuum limit of dedicated large-volume lattice studies, so once again the correlation of these errors are unknown.We conservatively take this correlation to lie somewhere between 0% and 100% and include this in the same two-point systematic as the other continuum limit correlations.
Taking into account both two-point systematics, we obtain four correlation matrices, each with an associated statistical uncertainty, By performing the fits described in the main text with each bootstrap sample for each of the four representative covariance matrices, we obtain the uncertainty in the rescaling factors and observed χ 2 values arising from our uncertainty in the lattice covariance matrix.
It is important to note that the statistical covariance matrices, C stat , used in the histogramming procedure, are determined from jackknifes over N J = 48 binned "measurements".Therefore, the statistical uncertainty on these covariance matrices is determined by N J , rather than by the number of independent samples used in the determination of the central values.This number is smaller than the total number of independent samples available in the dataset, even accounting for autocorrelations.This means that with the available data, a different sampling procedure (in particular, a larger number of bins) could produce an estimate of the statistical covariance with smaller statistical errors.However, this conservative approach to binning ensures that the systematic errors arising from autocorrelations within the data are negligible compared to these statistical errors.In addition, this choice provides consistency with the analysis of Ref. [6].As the results of Sec.V show, the statistical error this conservative binning produces is still sufficiently small that it does not quantitatively affect the conclusions drawn in the present work.

Appendix C: Weighted average procedure in presence of uncertainties on uncertainties and correlations
The minimum of the χ 2 function in Eq. ( 19), with respect, to γ is obtained for which represents a weighted mean of the input γj values.However, it is known that this approach can yield biased results (see e.g.Ref. [73]).
According to the Gauss-Markov theorem, this χ 2 minimization is equivalent to the minimization of the variance of a weighted average, with the sum of weights constrained to unity (see e.g.Refs.[97][98][99]).It is for this reason that this method is also called "BLUE", i.e.Best Linear Unbiased Estimate (see e.g.Ref. [97] and references therein).However, the word "unbiased" in the name above needs to be interpreted with care, as it implicitly involves several assumptions about the statistical model being employed.Indeed, in presence of strong correlations among the inputs, the weights in Eq. (C1) can be negative or larger than unity, yielding an average value that can be even outside the range of the input values, which may be problematic.Chapter 7 of Ref. [98] points out that the measurements are likely to lie on the same side of the true value in the presence of strong positive correlations.However, it is also indicated that this interpretation is subject to the implicit assumptions about the knowledge of the covariance matrix, as discussed in more detail below.
Ref. [73] points out that a bias in this procedure is caused by the input covariance matrix because of "the linearization on which the usual error propagation relies".While the χ 2 definition of Eq. ( 19) treats all the uncertainties as absolute (i.e.corresponding to additive effects), at least some of them should be treated as relative 12 (i.e.corresponding to multiplicative effects) [74].
Furthermore, another problem of this approach is that Eq. ( 19) involves an ill-behaved inversion of the covariance matrix.Indeed, a small change in the covariance matrix can have an important impact on its inverse, hence on the weights and the resulting average.Without special techniques, the inverse cannot be computed, nor can the χ 2 minimization be performed [74].While modern numerical algorithms allow to address the technical aspect of this matter, some more fundamental problems related to the evaluation of the covariance matrix still persist (see below).
In Ref. [73] it is discussed that a possible solution to these problems consist in introducing, in the χ 2 definition, one nuisance parameter for each correlated systematic uncertainty, together with the corresponding constraint terms.Doing so, for a normalization uncertainty, a scaling factor is applied to both data and uncorrelated uncertainties.An equivalent approach introduces the corresponding nuisance parameter as a scaling factor on the theoretical prediction [74,100].An alternative method for the treatment of relative uncertainties consists in scaling them to the fitted value minimizing the χ 2 , re-deriving the corresponding covariance matrix and iterating (see Refs. [73,101] and references therein) 13 .However we note that the rescaling of the uncertainties has little relevance in the current use-cases, when there's little difference between the γj values to be combined or when one of them is more precisely determined than the others and dominates the average (see Sec. V B).It is also worth recalling that, if all the nuisance parameters correspond to (are treated as) absolute systematic uncertainties, each of them acting as a coherent additive shift of the theoretical prediction, the approach using nuisance parameters is equivalent to the one of the "standard" χ 2 with correlations from Eq. ( 19) (see Refs. [102][103][104][105][106][107]).
The approach using nuisance parameters in the fit often leads to reduced systematic uncertainties, as they are constrained by the χ 2 minimization.However, this reduction assumes and strongly relies on a perfect knowledge of the phase-space dependence and correlations of the systematic uncertainties, an information explicitly used in the χ 2 definition.The same assumptions and features are implicit for the "standard" χ 2 with correlations in Eqs. ( 12), ( 19), (E8) and (E10), in particular for minimizing the variance of the average, as discussed above.Indeed, the need for precisely known covariance matrices, mentioned in Refs.[26,108], becomes explicit in the equations above.
Building upon remarks made in the context of ATLAS jet performance and cross-section studies [109][110][111] (see also the earlier discussion in Ref. [112]), these strong assumptions were also questioned for the e + e − annihilation data used for the evaluation of a LO-HVP µ , since there are clear indications that the amplitudes of systematic uncertainties and their correlations are generally impacted by uncertainties themselves [3,5,113]. 14Indeed, the size of the systematic uncertainties of the various measurements, the correlations of a systematic uncertainty component impacting several measurements, as well as the correlations between the systematic uncertainties impacting a given measurement are never really measured, but rather estimated.For example, the tensions between the input measurements observed in several channels (e.g. between BABAR and KLOE in the π + π − channel) are a direct indication of underestimated uncertainties for the measurements and motivated the inclusion of an extra systematic uncertainty for the resulting dispersive integrals [3] (see also Appendix A).Such two-point systematic uncertainty, although necessary when using the currently available inputs, has unknown correlations between different phase-space regions.Therefore, it comes at the price of bringing some fundamental limitations for the treatment of correlations in Eqs. ( 12), ( 19), (E8) and (E10).The recent CMD-3 measurement in the π + π − channel [13], which shows some significant tension with the previous precise measurements of the same channel, emphasizes even further the existence and the need for a careful treatment of the uncertainties on uncertainties.Furthermore, such limitations in the statistical treatment also originate from the uncertainties on uncertainties and on correlations for the lattice QCD calculations (see Appendix B).Therefore, the statistical procedures employing such inputs (or other quantities derived using them) should avoid overestimating the precision with which the uncertainties and their correlations are known, a conservative uncertainty treatment being preferable.This is especially important since the uncertainties on the uncertainties and on their correlations are not available in the current publications of hadronic spectra.
In order to address the limitations related to the uncertainties on the systematic uncertainties and on their correlations, we use a procedure which consists in performing a weighted average, with the weights of the input γj values being proportional to the inverse of their squared total uncertainties.This is equivalent to using only the diagonal elements of the covariance matrices (i.e.setting the other elements of the matrices to zero) in Eq. (C1) when deriving the weights.The weights used in this procedure are in the range between 0 and 1, and are on the larger end for more precise contributions.At the same time, the full information of the input uncertainties, with their correlations, is propagated to the result of the weighted average and is used in the χ 2 when assessing the level of agreement between the inputs and the average value.A similar fitting procedure, using only the diagonal uncertainties in the χ 2 followed by a propagation of the full set of uncertainties and their correlations, has been successfully used in Ref. [3] for combining the experimental measurements in the e + e − → π + π − channel. 15ppendix D: Stability of the results of Sec.V B with respect to the averaging procedure In this appendix we discuss the use of an alternative averaging procedure to the one employed in Sec.V B, to obtain the results for changes to the experimental Rratio that could explain the differences between the lattice and data-driven results for a LO-HVP µ , a LO-HVP µ,win and δ(∆ (5) had α).In that section, we perform the average defined in Eq. (C1).That is, when we derive the averaging weights, we replace the full covariance matrix, C γ lat + C γ R , obtained by a linear propagation of lattice and experimental R-ratio uncertainties and correlations, by its diagonal form with no correlation.Here we consider what the effect of including these correlations is on our results.
Generically, the use of the full covariance matrices for determining the averaging weights in Eq. (C1) changes these weights significantly, some becoming negative and others larger than unity (see Sec. IV B).However here, their deviations from zero or from unity are observed to be smaller than 0.06.Furthermore, in such cases either the γj values from Eq. ( 18) do not differ much 16 or one of them is more precise and dominates the average.Thus, the values of the average γ = γ 1 are similar for the two averaging methods, within 12% in the most extreme case, and in most cases within a few percent or less.
Even if this alternative approach minimizes the χ 2 in Eq. ( 19) (and the uncertainty of the rescaling percentage δ 1 , which is propagated from the covariance matrices of the lattice QCD and dispersive results), it only reduces it by 15% (9%) or less, compared to the values obtained with weights proportional to the inverse of the γj uncertainties squared, which are used in the nominal approach to produce Table II.Therefore, the p-values obtained with the two approaches are also similar, leading to similar conclusions about the level of compatibility between the lattice and dispersive moment integrals, for all the normalization shift scenarios considered in this paper.
The statistical variance and quantiles of the rescaling percentage δ 1 , obtained from the bootstrap variations of the lattice covariance matrices, are most often enhanced (by up to 540%) in the alternative averaging method compared to the one used to obtain the results given in the main text, while in some rare cases they are reduced (by up to 87%).At the same time, the changes of δ 1 under systematic variations of the lattice covariance matrix are also generally broader for the alternative averaging method.
Concerning the statistical variance and quantiles of the uncertainty on δ 1 (propagated from the covariance matrices of the lattice QCD and dispersive results), they are most often enhanced (by up to 133%) in the alternative averaging method compared to the nominal one, while in some rare cases they are reduced (by up to 99%).In all cases, as for the nominal method, the uncertainty of the rescaling percentage δ 1 remains precisely determined when using the alternative averaging method.
On the other hand, the variance and quantiles of the χ 2 values, due to the same bootstrap fluctuations as above, are generally reduced when using the alternative averaging method (by up to 44%), while in some rare scenarios they are slightly enhanced (by less than 3%).The range covered by the χ 2 values under systematic variations of the lattice covariance matrix is generally somewhat narrower for the alternative averaging method.Some of the features described above, for the comparison between the two averaging methods, are illustrated by comparisons between the plots in the first two rows of Figs.4-7.These features are consistent with the fact that while both procedures are sensitive to the uncertainty on the uncertainty in the determination of the weights, only the alternative approach, which uses the full covariance matrices for determining the averaging weights, relies on restricted √ s regions, below 0.63 GeV or above 1.8 or 3.0 GeV.In that case, the shape differences among the kernels of the various moment integrals have a more important impact for Eq. ( 18) and the differences between the γj values can be as large as 30 − 80%.
correlations for computing the weights to minimize the χ 2 and the uncertainty of the rescaling percentage δ 1 .
Appendix E: Going beyond the rescaling in a single √ s-interval To go beyond the simple scenarios and interpretations of Appendix C, more window observables are required from the lattice, with the corresponding ones from the experimental R-ratio, split up into √ s-intervals.These can be complemented by other quantities related to the HVP, such as the hadronic vacuum polarization function Π(Q 2 ) at various Q 2 accessible to the lattice.Again, all of these quantities must be obtained with a full description of the relevant covariances.

General formalism for comparing and possibly combining lattice and data-driven results
Consider a function M R (s, α i ), that depends on s and on parameters α i .It is chosen so as to provide a good description of the experimental R-ratio, or a re-binned version of it, once the α i are appropriately adjusted.Now, consider a modified model, M lat (s, α i , β j ), that depends on the additional parameters β j .This new model represents a guess at how the experimental R-ratio could be modified to reproduce the lattice observables of interest, once integrated over s with the appropriate weights.We then define the ratio of these two models as: which is chosen not to depend on the α i parameters. 17or instance, these functions can be physics-motivated models, spline-based parametrizations, orthogonal polynomials, histograms of the bins used for the experimental R-ratio or with some different binning, etc.For either the R-ratio or lattice model functions, we define the integrals of interest as where the α i and β j model parameters are kept implicit.
As above, I b can be either some restricted √ s-interval or it can cover the full range from threshold to ∞.In the latter case we drop the index b for the integral label, which also keeps the notations consistent with e.g.Eq. (15).Note that if M R/lat is a histogram, then the integral in Eq. (E2) becomes a sum over bins.
For the lattice constraints of the model integrals one can write the χ 2 function as: while, for the constraints from the experimental R-ratio, one can use approaches based on different χ 2 definitions, such as where, as above, we indicate whether the covariance matrix of the a j or a jb moment integrals correspond to those of the lattice (data-driven) approach via the subscript "lat" ("R").In addition, we call C R,bin the covariance matrix that corresponds to (a possibly rebinned version of) the experimental R-ratio.Indeed, Eq. (E4) concerns dispersive integrals on the full s range, Eq. (E5) includes further information on the shape of the R-ratio through the integrals computed on various I b intervals, while Eq.(E6) uses directly the experimental R-ratio (without any further convolution), either within its original bins or after merging the bins within larger intervals.Based on these χ 2 functions, one can consider several approaches for comparing and/or combining the R-ratio and lattice results.
A first approach can consist in constraining the parameters α i of M R (s, α i ) by minimizing one of the χ 2 R functions in Eqs.(E4)-(E6), followed by a determination of the β j parameters of M lat (s, α i , β j ) through a minimization of χ 2 lat from Eq. (E3).In this two-step approach, the α i parameters are held fixed when the second fit is performed and therefore they are not constrained by lattice results.Still, their uncertainties and correlations are propagated from the first step to the second, e.g. by replacing C lat in Eq. (E3) with the appropriate covariance matrix.
A second approach can consist in simultaneously constraining the α i and β j parameters of the model, by minimizing χ 2 lat + χ 2 R , with either of the definitions in Eqs.(E4)-(E6) for the contribution from the experimental R-ratio.We note that this minimum is generally lower than the χ 2 lat + χ 2 R value obtained for the models fitted in the two-step approach discussed above.While the minimization of the χ 2 defined in Eq. ( 12), according to Eq. ( 13) and Eq. ( 14), represents a combination of the lattice and R-ratio results in the space of the a j integrals, without any extra constraints, the approach of fitting the α i and β j parameters by minimizing χ 2 lat + χ 2 R represents a similar combination in the R-ratio space, under the constraint of the M lat (s, α i , β j ) and M R (s, α i ) models, respectively.
In a third approach, one can first minimize χ 2 lat from Eq. (E3) with respect to the free parameters of the model M lat (s, α i , β j ), obtaining αi and βj .One then considers a modified version of that lattice model which accounts for possible differences with the R-ratio.Using the notation of Eq. (E1), we define this model as Mlat (s, αi ) ≡ M lat (s, αi , βj )/F lat÷R (s, βj ).This model is subsequently binned via averaging within chosen √ sintervals.We call the resulting histogram M bin lat (s, αi ) and its covariance matrix Ĉlat,bin .In addition, we consider a model M R (s, α i ) for the R-ratio and its version MR (s, α i ), possibly re-binned to match the binning of M bin lat (s, αi ).The parameters α i are to be fitted simultaneously to the experimental R-ratio and M bin lat (s, αi ), by minimizing: where the first sum runs over the lattice bins and the second, over the ones of the experimental R-ratio.This effectively combines the lattice and experimental R-ratio results.
In this third approach, consider the case in which M lat is taken to be a binned histogram with the number of bins equal to the number of (linearly independent) lattice moment integrals, the count of each bin being a free parameter to be determined.Then the minimization of χ 2 lat of Eq. (E3) is equivalent to an "unregularized" unfolding towards the space of the R-ratio (in the sense that no other regularization besides the binning choice itself is being applied).Indeed, this is equivalent to solving a system of constraints through a matrix inversion approach (see e.g.Chapter 11 of Ref. [98]).However, the parametrization of M lat (s, α i , β j ) can be used to inject an effective regularization into the inverse problem.The unregularized approach avoids biases related to the matrix inversion.However, it may induce large variances and strong (anti-)correlations between the determined quantities (i.e.fitted R-ratio values in various bins) and hence a large hierarchy among the eigenvalues of the corresponding covariance matrix.In that case, any further use of such results would require a very precise determination of the covariance matrix which, in turn, would necessitate very precise knowledge of the input covariances.Alternatively, it is generally preferable to use regularized unfolding methods, aiming for a trade off between moderate statistical variances and systematic biases induced through the regularization. 18 It is to be noted that the three approaches described above, with their possible variants, imply different hypotheses about the validity of the fitted models and the way in which their parameters are constrained.The resulting χ 2 /ndof values represent tests of these hypotheses, providing information about the compatibility between these models and the input data, as well as about the compatibility between the lattice and R-ratio inputs themselves.In the presence of tensions, one can consider enhancing the uncertainties of the combination result, following e.g. the approaches from Refs.[3,[70][71][72], summarized in Appendix A. This can be achieved by using the information from some (partial/local) χ 2 and/or adding some extra uncertainties to account for systematic deviations between the inputs.
Moreover, while Eqs.(E3)-(E6) are written using the full covariance matrices, the discussion of Appendix C, concerning the treatment of the uncertainties on the uncertainties and on the correlations is relevant here too and has to be implemented accordingly.For instance, one can perform the χ 2 fits for determining the model parameters using only the diagonal elements of the covariance matrices, while the full information on the uncertainties and their correlations is propagated to evaluate the uncertainties of the model parameters, their correlations and the fit quality.To make contact with what was done in the main body of the paper, we consider here the rescaling in a single √ s-interval that was introduced in Sec.IV B and employed in Sec.V B. It can be retrieved as a special case of the general approaches discussed above.Indeed, one can consider a "model" M R (s, α i ) that consists in a histogram with the same binning as the experimental Rratio itself.The α i are parameters to be determined and correspond to the bin counts.Employing the first, twostep approach described below Eq. (E6), with χ 2 R given by χ 2 R,3 from that equation, the first minimization step yields the bin counts with χ 2 R = 0. Moreover it propagates the covariance matrix of the experimental R-ratio to the histogram of M R (s, α i ), i.e. to the parameters α i .The model is completed with the function F lat÷R (s, γ), which is equal to a rescaling factor γ in a subset of the √ sintervals, while it is equal to unity for the complementary subset of bins.This yields a model M lat (s, α i , γ) that implements the prescription from Eq. (17).In particular, this model ensures that the ratios of integrals a MR jb /a MR kb are to be determined, allowing us to also address the worries that are sometimes expressed concerning the ill-defined nature of the problem when using χ 2 -based approaches (see e.g.Refs.[34-36, 38, 41, 42]).
and a M lat jb /a M lat kb are equal to a R jb /a R kb .In turn, this guarantees that these ratios are independent of the parameter γ, for any j and k integrals and any interval I b .In other words, these ratios of integrals are driven by the corresponding kernels convoluted with (i.e.averaged over) the shape of the experimental R-ratio, while the changes of normalization through the factor γ impact coherently all the a M lat jb integrals on the rescaled √ s-interval.Eq. (E3), with the model M lat (s, α i , γ) and C lat completed with a propagation of the covariance matrix of the parameters α i , yields Eq. ( 19) (after a change of variable similar to Eq. ( 18), together with the corresponding linear uncertainty propagation).
The approach discussed above can be generalized by applying the constraint from Eq. ( 16) at the level of the corresponding R-ratio and lattice models, with multiple rescaling factors γ b applied in the √ s-intervals I b .This corresponds to minimizing provided that the number of linearly independent moment integrals is greater or equal to the number of fitted γ b rescaling factors.In this equation, C lat + CR is the covariance matrix corresponding to a lat j − b γ b a R jb , depending hence also on the γ b parameters.As above, this implies assigning a pre-fit uncertainty to the M R (s, α i ) model, which is then propagated via the corresponding covariance matrix. 19In addition, if the I b intervals correspond to the original experimental R-ratio bins and their number is equal to that of linearly independent lattice observables, a lat j , then the rescaling hypothesis becomes exact.In that case, one can have as many normalization rescaling factors as there are bins, allowing for a complete, unregularized reconstruction of the lattice R-ratio with the same granularity as that of the experimental R-ratio.Please see the discussion after Eq. (E7) for more details.It is also worth noting that in the case when all the γ b = 1, the χ 2 from Eq. (E8) matches the one of Eq. ( 14), providing a direct compatibility check between the lattice and the dispersive results. 19These methods can also be seen as fits of normalization factors for templates of moment integrals derived through the dispersive approach, used to describe the corresponding lattice results, as well as possible.For a discussion on the template fitting methods, with the various approaches for treating the uncertainties of the fitted distributions and of the templates themselves, see e.g.

Single or multiple rescaling factors for fitted moment integrals
An alternative to Eq. (E8), which includes multiple rescaling factors, γ b , is given by the model that is assumed to describe the inputs, a lat j and a R jb , within the uncertainties.The second line provides a trivial model for the R-ratio contributions, a R jb , from the √ s-intervals, I b , defined after Eq. ( 15) and the first a rescaling model for the corresponding lattice observables defined in Sec.III.Of course, the number of lattice observables, a lat j , must be larger than the number of free rescaling parameters.This is a special case of Eq. (E1), as given in Eqs.(E3) and (E5) without, however, the constraints of the models, M R/lat , for the lattice and R-ratio spectral functions.
Thus, to determine the parameters of this model, taking into account all correlations, one can minimize with respect to the γ b and a jb .Again, lattice results are assumed to have no correlations with those obtained from the experimental R-ratio.The uncertainties on the parameters γ b and a jb can be obtained from the Hessian or by using pseudo-experiments.One can also consider the p-value obtained from the minimum value of the χ 2 and the system's number of degrees of freedom, which is the number of lattice window observables minus the number of rescaling parameters.This p-value indicates the consistency of the available lattice and R-ratio input with the rescaling hypothesis above.Now, suppose that the rescaling model of Eq. (E9) is well justified theoretically and that the p-value of the fit is acceptable.Then, the resulting parameter values provide a description of the lattice and R-ratio observables, via the right-hand sides of Eq. (E9), that combines the information contained in the a lat i and a R jb .When the number of observables and intervals increases, so will (anti-)correlations between them, leading to large variances on the results, well known in inverse problems.Of course, if the number of rescaling parameters is the same as the number of lattice observables, then the system of equations of Eq. (E9) can be solved exactly, and one obtains a jb = a R jb as well as the values of the γ b obtained by minimizing Eq. (E8), in the same situation.
Beyond the possibility of combining lattice and R-ratio results for a lat j and a R jb , via the model of Eq. (E9), the approaches of Sec.IV B and of Eqs.(E9) and (E10) differ in that the former is based on a model that is linear in the fitted parameters (Eq.( 18)), with a linear propagation of uncertainties, while the latter refers to a nonlinear model (Eq.(E9)), but with no linear propagation of uncertainties.They are expected to be equivalent in the limit that the rescaling factors are such that the |γ b − 1| are small and as long as the relative uncertainties on the a R jb are also small.We have checked this analytically and numerically, in the case of a single rescaling factor γ (either fitted or fixed to 1), for the results presented in Sec.V.The non-linear effects in the uncertainty propagation mentioned above were also studied and found to be numerically small.The non-linearities in the deviation from 1 of the rescaling parameter, γ, are also generally small, somewhat larger when the model does not describe the lattice and R-ratio results well.
The two approaches also answer slightly different questions.The one of Appendix E 2 corresponds to a rescaling of the dispersive integrals derived from the experimental R-ratio in different √ s sub-intervals I b , while the one presented here, of a rescaling of the fitted a jb , which may deviate from the input a R jb .
Moreover, the approach described here provides an unregularized reconstruction of both the lattice observables, via the first line of Eq. (E9), and of the timelike contributions of √ s-intervals, I b , as given in the second line of Eq. (E9).The resolution with which this sort of reconstruction can be performed is usually limited by (anti-)correlations between the various observables that generally leads to large variances on the results.Note that, when all the γ b factors are set to unity, one can perform the change of variable a j b = a j − b̸ = b a jb in Eq. (E10), for an arbitrary index b.The covariance matrix CR , resulting from the propagation of C R through this change of variable, can be decomposed into a set of uncertainties that are either uncorrelated (σ j(b) ) or that correlate (s with the index l running over the corresponding independent uncertainty components [107].This enables an equivalent re-writing of the χ 2 function from Eq. (E10) 20 In quantities such as a R j(b) , the parentheses indicate that we are considering both a R j and a R jb with b ̸ = b. as using nuisance parameters, s l j(b) , to account for the correlated uncertainties [102][103][104][105][106][107]. 21While the a j parameters are constrained by both the lattice and the dispersive inputs, the a j,b̸ = b parameters remaining after the change of variable only enter the dispersive part of the χ 2 function.The minimisation of this χ 2 function with respect to a j,b̸ = b allows us to exactly cancel the corresponding contribution to the χ 2 , regardless of the contributions of the shifts of the nuisance parameters away from zero, which are parameterized by the β l and induced by the a R j and a lat j inputs.Switching back from the representation in terms of nuisance parameters to the one with the covariance matrix, this time in the subspace of a j parameters, allows to show that the minimum of Eq. (E10) with respect to a j,b̸ = b reduces to Eq. ( 12), while their global minimum is provided by Eq. ( 14).This also implies that for γ b = 1 the minimum of the χ 2 function from Eq. (E10) is stable with respect to the choice of the I b intervals, a feature also observed numerically.

Parameterizing possible shape changes of the R-ratio suggested by lattice results
Once sufficient lattice input is available, one can refine the analysis in other ways.Since our results suggest that the ρ peak could be responsible for the disagreement with the data-driven approach, one could imagine focusing on that region and subdividing it into smaller √ s-intervals.Then, with the methods described in Appendices E 2 and E 3, one could reconstruct, via rescaling parameters, the modifications to the experimental R-ratio suggested by the lattice observables.

a. Using complete sets of polynomials
Because we are interested in understanding how a tension between lattice and data-driven results could suggest a necessary alteration of the normalization, slope, 21 In Eq. (E12), the parentheses around the b and b indices indicate that the sum in the second term of the equation runs over the a R j and a R jb .
curvature, etc. of the R-ratio in a given √ s-interval, it also makes sense to parameterize the deviation from the experimental R-ratio as a sum of Legendre polynomials, P ℓ , that form a complete basis on that interval.Thus for simplicity, we choose the model M R (s, α i ) to correspond to the histogram of the experimental R-ratio, with the original binning and the α i parameters tuned to the corresponding bin counts.We also choose the function F lat÷R of Eq. (E1) to have the following form:22 2 ] and I 2 is the complementary interval on the real axis.
Then we write: where a M lat jb , b = 1, 2, are defined via Eq.(E2), with the appropriate kernel Kj (s).The parameters of the model are determined via the first approach described in Appendix E 1.The coefficients, β ℓ of the Legendre polynomials are then determined by minimizing χ 2 lat of Eq. (E3).The resulting model, M lat (s, α i , β ℓ ), serves as a regularization of the unfolding of the lattice R-ratio, in the sense discussed at the end of Appendix E 1, provided that the number of fitted Legendre coefficients is smaller than the number of moment integrals being considered.
Note that the rescaling hypothesis of Sec.IV B is a special case of the one considered here, with ℓ restricted to the value 0 and δ = β 0 − 1.

b. Using physics-driven models
Up until this point, none of the considered modifications to the experimental R-ratio incorporate theoretical constraints which the spectral function is known to obey.Parametrizations that do can be found in Refs.[3,24,[119][120][121][122][123].
Here, for concreteness, we focus on the parametrization given in Ref. [24].It relates the R-ratio to the π + electromagnetic form factor in the timelike region and it is the latter which is actually parametrized.This parametrization satisfies many requirements of analyticity and unitarity.It also includes 2 and 3-pion channels, as well as inelastic corrections via a conformal polynomial, with a threshold dictated by phenomenology.The authors estimate that this parametrization is valid from the two-pion threshold to √ s ≃ 1 GeV.To fit it to the measured Rratio, the latter must be corrected for final-state radiation [24].
Because we do not actually analyze results with this parametrization here, we do not describe it in detail.Its expression is given in Ref. [24].However, we make a few comments on its parameters and how it may be used to describe a modification to the measured R-ratio compatible with lattice results.
The first set of parameters is the two-pion phase shift at √ s 0 = 0.8 GeV and √ s 1 = 1.15 GeV.The second set of parameters describes ρ-ω mixing, in terms of the mass and width of the ω, as well as of a mixing parameter.
The last set of parameters describes inelasticity above √ s in = M ω + M π 0 , with typically 4 parameters.That constitutes a set of nine parameters that we label p k .In Ref. [25], the authors perform a study of the modifications to the measured R-ratio implied by a value of a LO-HVP µ that differs from its data-driven prediction.However, the assumptions behind that study, and the questions which it answers are quite different from those of the approach proposed here.
In Ref. [25] the authors first perform a fit of the parametrization of Ref. [24] to all e + e − → π + π − measurements available at the time, fixing its nine parameters.This allows them to predict the two-pion contribution to a number of quantities, including a LO-HVP µ , in the data-driven approach.Then they release some of the nine parameters, keeping the others fixed to their data-driven values, while performing the fit once more, with the constraint that the predicted value of a ππ µ takes on a precise, prescribed value.This value is allowed to vary from its nominal, data-driven value to the value of a ππ µ that one obtains by assuming that the full difference between a lattice result for a LO-HVP µ , close to that of Ref. [6], and its data-driven counterpart can be ascribed to the two-pion contribution.In particular, the approach allows the authors to determine the ∆χ 2 between a scenario in which a ππ µ takes on its "lattice" versus its nominal value.In turn, the value of the ∆χ 2 indicates how compatible different values of a ππ µ are with e + e − → π + π − data and known constraints on the pion electromagnetic form factor.The approach also allows the authors to present correlations between a ππ µ values and predictions for the two-pion contribution to other observables, such as the pion electromagnetic radius or the running of α.
Here, in a sense that we now describe, we address the problem the other way around.As discussed throughout this article, we choose to consider the measured R-ratio, or quantities derived from it, and as many observables as possible that are computed in lattice QCD and that are related to HVP.The latter include, for instance, the contributions to a LO-HVP µ in various time windows, the running of α in certain spacelike intervals, etc. Concerning the spectral function aspect of the comparison, we not only consider the √ s-interval I 1 = [2M π ± , 1 GeV], but also its complementary interval I 2 .As for all the comparisons between the two approaches that we have looked upon until now, knowledge of the correlations between all quantities considered must be determined as precisely as possible.Then, in the notation of Appendix E 1, we call M R (s, α i ) the model for the R-ratio.To be more specific, in I 1 one can consider M R (s, α i ) to be the model of Ref. [24], so that nine of the α i correspond to the parameters p k discussed above.We continue calling α i the other parameters of the model.To describe the R-ratio in I 1 , only the p k are needed.The α i are involved in describing the R-ratio in I 2 .In that interval, M R (s, {p k , α i }) can be chosen to be any of the models discussed above for quantities related to the R-ratio, e.g. the histogram of the measured spectral function, with parameters α i corresponding to the counts in the original, measurement bins.
Now we consider a model for the candidate lattice Rratio, M lat (s, α i , β j ).Instead of defining it via Eq.(E1), we take it to be the same model as M R (s, {p k , α i }) in I 1 , but where some of the parameters β j correspond to possible shifts, δp l , of a subset of the parameters p k of the model of Ref. [24].Thus, as for M R (s, {p k , α i }), we rewrite the lattice R-ratio model M lat (s, {p k , α i }, {δp l , β j }), where the parameters α i are those needed to describe the measured R-ratio in I 2 , via M R , and the β j , those that are used to express the differences between the measured and lattice spectral functions in that same interval.For instance, the β j could correspond to rescaling factors of the measured R-ratio in sub-intervals of I 2 , discussed in Appendix E, or Legendre polynomial modifications, as discussed in Appendix E 4.
The next step consists in following the second approach described in Appendix E 1, in which the parameters p k , δp l , α i and β j are simultaneously determined by minimizing χ 2 lat + χ 2 R , where χ 2 lat is given by Eq. (E3) and χ 2 R is one of the χ 2 given in Eqs.(E4) and (E5).In this physics-motivated approach, both the measured and the candidate lattice R-ratios are described by a model that incorporates as many physical constraints as desired in interval I 1 .Moreover, the possible differences between the data-driven and lattice HVP are expressed in terms of parameters which have a physical meaning.For instance, suppose that the δp l only correspond to modifications in the two-pion phase shift at √ s 0 = 0.8 GeV and √ s 1 = 1.15 GeV.If the minimization of χ 2 lat + χ 2 R gives a good χ 2 value and at least one of the δp l is not consistent with zero, then we can argue that, in the interval I 1 , the lattice and data-driven approaches differ in their prediction for the phase of the pion form factor.In addition, if one has reasons to believe that M lat correctly describes the differences between the lattice and data-driven approaches and one knows which of the two approaches correctly describes HVP, then the corresponding model can be viewed as a combination of the data-driven and lattice results, and can be used to make improved predictions for quantities that involve HVP.
We leave the comparison of lattice and data-driven results via physics-motivated models, as described here, for future work.In Ref. [26] it was stated that several window integrals, with possibly narrower [t min , t max ] ranges compared to the ones studied up to now, could be employed to perform a more detailed comparison of the dispersive and lattice approaches.While such a larger set of window moments generically enhances the amount of available information, this increase will eventually be limited by the (anti-)correlations among the windows, for both the dispersive and lattice determinations.
In the following studies we use a "blinding" approach, in the sense of only communicating the information on the dispersive uncertainties and their correlations for the moment integrals, but not the corresponding nominal values.This is motivated by the fact that these moments have not yet been determined via lattice QCD and we wish to preserve the possibility of doing so without any risk of bias.This is different to the approach followed in Ref. [26] that quotes nominal values for moment integrals derived through the dispersive approach (obtained with different data combination methodologies from those used here), even when these integrals have not yet been computed via lattice QCD.We believe that going forward, all groups working on this problem should adopt our more conservative approach.Employing "blinding" approaches is very important when working at the level of precision involved in these studies.This is especially the case in the current context of the existing tensions between the dispersive approaches and the lattice calculations, as well as among the various experimental results themselves.
We consider first the seven window moments discussed in Ref. [26], then a set of nine window moments where the [1.6, ∞[ fm window is split into three, with extra thresholds at 2.6 and 4 fm, and finally a set of nineteen moments including ten extra ∆α (5) had (q 2 ) results, with spacelike q 2 varying from −10 GeV 2 to −1 GeV 2 .Table V shows the correlation matrix of these nineteen moment integrals, computed using the dispersive approach.It represents a first way of quantifying the amount of independent information available in these moments.Indeed, one can notice that moment integrals with similar kernel shapes (e.g.∆α (5) had (q 2 ) moments computed for similar q 2 values) are strongly correlated, while other moments have weaker correlations.This reflects the fact that differing kernel shapes supply a more effective way of exploiting the information provided by the timelike spectra.
A way of further quantifying the amount of independent information in the moment integrals is to evaluate the eigenvalues of the corresponding covariance matrix, since strong correlations typically imply the presence of small eigenvalues.Tables VI, VII and VIII show the dispersive uncertainties for the three sets of moments, as well as the eigenvalues of the corresponding covariance matrices, which typically show a fast drop over several orders of magnitude.It is worth noting that the uncertainties of the ∆α (5) had (q 2 ) and a LO-HVP µ,win integrals typically differ by several orders of magnitude, just because of the natural normalization of these moments, which is also reflected in the range covered by the corresponding eigenvalues of the covariance matrix.This effect, also present to some extent when considering separate series of a LO-HVP µ,win and ∆α (5) had (q 2 ) integrals, interferes with the intent of quantifying the amount of independent information in the moments.
To take this into account, we also consider the eigenvalues of the correlation matrices of the moment integrals, as well as of the covariance matrices normalized by the nominal values of the moments (i.e.[C R ] ij / a R i • a R j ).These are also displayed in Tables VI, VII and VIII 23 24 .For the set of seven a LO-HVP µ,win dispersive moment integrals in Table VI, these eigenvalues cover seven orders of magnitude, which reflects the strong correlations among some of these moments.For the three sets of moments, while the largest eigenvalues are enhanced by a factor of 2 − 3, as more moments are added, the smallest eigenvalues drop by more and more orders of magnitude.This indicates that the amount of independent information (i.e. of independent degrees of freedom) is less than the number of added extra moments.Indeed, when comparing the ratios of the various eigenvalues with the largest one, a reduction of the eigenvalues by two or four orders of magnitude for Table VI allows us to accommodate typically one extra moment in Table VII.At the same time, when including the extra ∆α (5) had (q 2 ) moments in Table VIII, the potential gain in extra degrees of freedom seems rather marginal.This is related to the fact that the ∆α (5) had (q 2 ) moments considered here are rather strongly correlated among themselves and have rather strong correlations with at least one of the a LO-HVP µ,win moments (see Table V).One can attempt to provide more quantitative conclusions concerning the number of independent degrees of freedom that are available, in the data-driven approach, for a comparison with lattice results.For this purpose, one can consider the eigenvalues of the correlation or normalized-covariance matrix, rescaled by the largest eigenvalue.The number of independent degrees of freedom can then be obtained by counting the number of such ratios above a chosen lower bound.This procedure is consistent with the idea that there is a limit to the precision with which we can obtain the uncertainty of moments of the measured R-ratio and linear combinations thereof.Such a limit approximately corresponds to that lower bound.TABLE V. Correlation matrix of ten ∆α (5) had (q 2 ) and nine a LO-HVP µ,win dispersive moment integrals computed for various q 2 values and [tmin, tmax] intervals with ∆ = 0.15 fm, respectively.

Moment integral
Correlation coefficients ∆α (5) had (−10 GeV 2 ) 1 ∆α (5) had (−9 GeV 2 ) 0.999 1 ∆α (5) had (−8 GeV 2 ) 0.999 0.999 1 ∆α        Here we choose this lower bound to be 10 −6 .With this bound we find that the seven window observables listed in Table VI contain a significant amount of linearly independent information.This number rises to eight when two more long-distance window observables are considered in Table VII.However, when ten values of ∆α (5) had (q 2 ) are added, in Table V, to the nine window observables of Table VII, no additional independent degrees of freedom appear.Thus, assuming that a lower bound of around 10 −6 is appropriate, it is reasonable to conclude that the number of linearly independent moments which can be obtained from the experimental R-ratio is less than ten for the types of observables considered here, which are the ones that should be computable, in the lattice approach, with sub-percent precision.
All these and other results presented in the paper point to the importance of determining the corresonding matrices in the lattice approach and performing studies similar to those presented here, when more moments become available in that approach.
had (q 2 ) and nine a LO-HVP µ,win dispersive moment integrals computed for various q 2 values and [tmin, tmax] intervals with ∆ = 0.15 fm, respectively.The eigenvalues of the corresponding covariance, correlation and normalized covariance matrices are also shown.

Eigenvalues Moment integral
Total uncertainty Covariance Correlation Normalized covariance ∆α

2 Z
) are observed for modifications of the R-ratio in the high-√ s intervals.However, while such modifications may have reasonable p-values with constraints from only a LO-HVP µ and a LO-HVP µ,win

µ
and a LO-HVP µ,win , is unfeasible.Instead, we prepared three one-dimensional histograms for the distributions of X = a LO-HVP µ , Y = a LO-HVP µ,win and Z = a LO-HVP µ + a LO-HVP µ,win B28) for C stat and C syst .Here, C stat and C syst are the statistical and systematic covariance matrices of a LO-HVP µ and a LO-HVP µ,win

FIG. 8 .
FIG. 8. Histograms for the estimation of the total, statistical and systematic lattice covariance matrices for the connected ud (top), the connected s (middle) and disconnected (bottom) contributions to a LO-HVP µ

TABLE II .
Results for various rescaling scenarios, using two observables (aLO-HVP µ and a LO-HVP µ,win

TABLE III .
Starting time for the taste corrections in the 12 systematic combinations used to estimate the correlation introduced by the taste correction.The first column contains the number of the systematic combination, the second column contains the starting time for the taste correction in a LO-HVP µ and the last column contains the starting time for the taste correction for aLO-HVP

TABLE IV .
Upper and lower bounds of the histograms used for determining the lattice covariance matrices in different flavor channels.