The Euclidean Adler Function and its Interplay with $\Delta\alpha^{\mathrm{had}}_{\mathrm{QED}}$ and $\alpha_s$

Three different approaches to precisely describe the Adler function in the Euclidean regime at around $2\, \mathrm{GeVs}$ are available: dispersion relations based on the hadronic production data in $e^+e^-$ annihilation, lattice simulations and perturbative QCD (pQCD). We make a comprehensive study of the perturbative approach, supplemented with the leading power corrections in the operator product expansion. All known contributions are included, with a careful assessment of uncertainties. The pQCD predictions are compared with the Adler functions extracted from $\Delta\alpha^{\mathrm{had}}_{\mathrm{QED}}(Q^2)$, using both the DHMZ compilation of $e^+e^-$ data and published lattice results. Taking as input the FLAG value of $\alpha_s$, the pQCD Adler function turns out to be in good agreement with the lattice data, while the dispersive results lie systematically below them. Finally, we explore the sensitivity to $\alpha_s$ of the direct comparison between the data-driven, lattice and QCD Euclidean Adler functions. The precision with which the renormalisation group equation can be tested is also evaluated.

1 Introduction Two-point functions are among the most basic objects that one can define within a Quantum Field Theory (QFT). Every possible action encodes associated distributions for them. From the phenomenological point of view, a particularly interesting one is the two-point correlation function of two vector neutral quark currents, Π µν ij (q) ≡ i d 4 x e −iqx 0|T (q i (x)γ µ q i (x)q j (0)γ ν q j (0)) |0 = (q µ q ν − g µν q 2 ) Π ij (s ≡ q 2 ) .
(1.1) Even in the partonic approximation, this is in the absence of quantum corrections and neglecting quark masses, these two-point correlation functions are divergent and depend on an arbitrary subtraction prescription. This approximation makes sense as the leading term of an Operator Product Expansion (OPE) [1][2][3] which is well defined for large Euclidean momenta, Q 2 = −q 2 and leads to 1 2) where N C is the number of quark colours. In order to avoid unphysical subtraction ambiguities, it is convenient to define the object of study in this work, the Euclidean Adler function [4]: which gives in the partonic limit The most important corrections in the Standard Model (SM) come from strong interactions. At large Euclidean momenta, far enough from the non-analytic behaviour in the distributions induced by hadrons, deviations from asymptotic freedom are described by perturbative QCD (pQCD), so that the leading corrections are simply given by δD ij (Q 2 ) = N C δ ij αs(Q 2 ) π . For light quarks i, j ≤ 3, D L ij (Q 2 ), the perturbative QCD description breaks down in the infrared region: the low-energy hadron dynamics is not well described by approximately-free quarks and gluons. Our knowledge about it at very low energies (there are no physical singularities at Q 2 = 0) implies that D ij (Q 2 → 0) = 0. 2 1 In the rest of the complex plane Πij(Q 2 ) can be obtained by analytic continuation. We take Arg(Q 2 ) ∈ [−π, π) and q 2 = e iπ Q 2 , so that Arg(q 2 ) ∈ [0, 2π). In this way, ImΠ(|Q 2 |e −iπ ) = 1 2i [Π(|Q 2 |e −iπ ) − Π(|Q 2 |e iπ )]. 2 Let us note that if one defines α eff s (Q 2 ) ≡ π D ii (Q 2 ) N C − 1 , then α eff s (0) = −π, which per-se does not tell us any new fundamental knowledge about strong interactions. Low-energy Effective Field Theories (EFTs), such as Chiral Perturbation Theory [5][6][7], give some nontrivial information about the infrared behaviour of D L ij (Q 2 ), but their predictive power is limited, especially at intermediate energies. For massive quarks, i, j > 3, D H ij (Q 2 ), a perturbative QCD description is known to give a precise description of the Adler function even in the neighbourhood of Q 2 → 0 [8][9][10][11], since quark masses regularize, at least up to a certain extent, the gluon singularities associated to infrared propagators.
Other powerful nonperturbative methods can also be used to obtain D ij (Q 2 ). Numerical simulations in a discretized space-time lattice [12] allow for a precise computation of the two-point functions at Euclidean momenta without relying on perturbation theory. Indeed huge efforts are recently being made to compute the electromagnetic correlator, where Q i is the electromagnetic charge of the associated quark in units of e (e.g. Q 1 = 2 3 ), since it plays a fundamental role in our understanding of the anomalous magnetic moment of the muon and in the so-called hadronic running of the QED coupling [13][14][15][16][17][18][19][20][21][22][23]. Currently the predictive power of lattice methods becomes severely limited as one goes above Q ∼ 2 − 3 GeVs due to discretization effects, leading to an interesting complementarity with respect to pQCD.
Similar motivations have increased the knowledge of Π(Q 2 ) obtained from another powerful nonperturbative method, the dispersive data-driven approach, which mainly uses electron-positron data to determine Π(Q 2 ) [24][25][26][27][28][29]. A well-known limitation in the current precision comes from a series of tensions involving electron-positron data. Besides the long-established discrepancy between e + e − → π + π − and τ − → π − π 0 ν τ data (invoking an isospin rotation) [30], there are significant tensions among different e + e − data sets (mainly KLOE vs BABAR, for the same 2π channel). Additionally, there is a clear discrepancy between the experimental value of (g−2) µ and the theoretical SM prediction obtained when using e + e − data to evaluate the Hadronic Vacuum Polarization (HVP) contribution (even after inflating uncertainties to account for the KLOE-BABAR tension). Finally, further tensions arise between e + e − data and lattice evaluations of both the hadronic running of the QED coupling and again the HVP contribution to (g − 2) µ [13,14,[31][32][33][34].
In this work we study analytically the Euclidean Adler function, in the Q 2 region where perturbation theory is expected to be valid, with the aim of comparing it to both the dispersive Adler function D(Q 2 ) obtained with the DHMZ compilation of data [28] and the one obtained from recently published lattice results for the hadronic running of the QED coupling [14].
On the one hand, assuming the validity of pQCD at a certain Euclidean momentum Q, one can check whether the description is consistent with the other approaches. On the other hand, assuming that the other approaches are correct, we have a uniquely clean window to learn about the onset of the asymptotic regime and the value of the associated QCD coupling. We start by introducing the overall theoretical framework to connect the different descriptions of the Euclidean Adler functions and the HVP. This is done in Sec. 2. In Sec. 3 we perform a comprehensive study of the perturbative Adler function in the regime we are interested in, combining many existing results. Then, in Sec. 4, we explain in detail how the data-driven and the lattice Adler functions are obtained, and in Sec. 5 we compare them with the perturbative expression. Power corrections are discussed in Sec. 6. Finally an exploration to the sensitivity of the comparison to the strong coupling and discussion about some possible fitting strategies can be found in Sec. 7. Conclusions and final remarks are presented in Sec. 8.

Theoretical framework
Let us start by connecting the needed descriptions and observables related to the HVP. The hadronic running of the QED coupling can be defined in terms of the electromagnetic correlator Π(Q 2 ) as follows, with α = α(0) = 1/137.035 999 084 (21) [35]. Since both the hadronic running of α(Q 2 ) and the Adler function are defined in terms of the electromagnetic correlator (cf. Eq. (1.3)), it is straightforward to relate both of them, Alternatively, On the other hand, it can be shown that the ratio of hadronic and muonic production cross sections in e + e − annihilation is directly related to the imaginary part of Π(Q 2 ). Formally, it is defined as R(s) ≡ 3s 4πα σ 0 (e + e − → hadrons (+γ)) = 12π ImΠ(Q 2 = s e −iπ ) , (2.4) where σ 0 (e + e − → hadrons (+γ)) refers to the so-called "bare" hadronic cross section in electron-positron annihilation, subtracting the vacuum polarization contribution to the photon propagator and the Initial State Radiation (and ISR-FSR interference) but including all Final State Radiation. Thus, in order to relate R(s) with D(Q 2 ) and ∆α had (Q 2 ), we have to relate the imaginary part of the analytic continuation of Π(Q 2 ) with Π(Q 2 ) itself. This is accomplished by using dispersion relations, which combine our knowledge of the analytic structure of the electromagnetic correlator, Π(Q 2 ), with its known asymptotic behaviour. Indeed, the partonic picture, see Eq. (1.2), is valid at very high energies and Π(Q 2 ) is an analytic function in the whole complex plane except for a branch cut starting at the hadronic threshold, 3 As a consequence, we can integrate Π(Q 2 ) in the complex plane along any contour as long as we avoid this cut. We may also weight the integral with some function , (2.5) which is analytic in the complex plane except for the two poles. Then, integrating along the contour in Fig. 1 and using the Residue Theorem, we get (2.6) We may separate the different contributions to the contour integral as (2.7) For |Q 2 0 | → ∞ one can easily show that the first term in the first line goes to zero by using the partonic description. Performing a change of variables and using the Schwarz reflection principle in the second term, one arrives at Starting from this expression, one can obtain dispersion relations for different objects by choosing different values of Q 2 1 and Q 2 2 . Choosing either Q 2 1 or Q 2 2 to be 0 leads to a dispersion relation for the correlator Π(Q 2 ), It is then straightforward to write an equation for the data-driven determination of the Euclidean Adler function using Eq. (2.10) and Eq. (2.4), (2.11) Finally, ∆α had (Q 2 ) is also related to R(s) through the dispersion relation of Eq. (2.9) and Eq. (2.4), . (2.12)

The perturbative Adler function
The perturbative Adler function has been the subject of many studies in different energy regimes, since it is directly linked to different precisely known inclusive observables [36]. At order α 2 s and beyond, the coefficients in the expansion depend on the renormalization scheme. For practical purposes, the MS scheme is usually adopted due to its computational simplicity. In the limit of n f massless quarks and no massive ones, the corresponding Adler functions are known up to (and including) five loops (i.e. order α 4 s ). The real world contains six quark flavours with a striking hierarchy of quark masses: where Λ QCD is the QCD scale, so that perturbative QCD does not make sense below it. In mass-independent renormalization schemes such as the MS, the perturbative series with six quark flavours does not give a very accurate approximation to the QCD Adler function at |Q 2 | m 2 t . This is a consequence of the lack of decoupling associated with this type of schemes, which leads to log (m 2 t /Q 2 ) factors that slow down and eventually break the perturbative series (e.g. see [37] for a pedagogical introduction to the problem). If one wants to keep track of the full quark-mass dependence of the Adler function [38], useful for example for the perturbative running of ∆α had (Q 2 ) from Q ∼ 2.5 GeV to Q = M Z [39], one possibility is employing a renormalization scheme that automatically performs the decoupling of heavy masses, such as MOM, at the cost of more complex calculations and less known perturbative corrections. Alternatively, one can still work in the MS scheme by introducing a series of QCD effective field theories with different number n f of massless quark flavours, which need to be matched at the corresponding quark-mass thresholds. The decoupling of heavy masses is then implemented by hand and the massless Adler function can be supplemented with the corresponding power-suppressed corrections from the heavy quark masses. The small contributions from the non-zero light-quark masses can be taken into account through perturbative expansions in powers of m 2 q /Q 2 . Since the other methods analyzed in this work are also more powerful below |Q 2 | ∼ 5 GeV 2 , we will focus on the Adler function D(Q 2 ) at |Q 2 | < 4m 2 c . In the next subsections we study the different perturbative Adler functions D ij (Q 2 ) separately, first focusing on those with i, j < 4 (light-quark contributions) and then on the heavy ones, i, j > 3 (heavyquark contributions). Eventually we add the leading QED corrections, put them together and also explain why the mixed (i < 4, j > 3 and vice-versa) contributions to D(Q 2 ) are very suppressed.

Leading massless contributions
For light flavours i, j < 4, the MS scheme with n f = 3 massless quarks gives an accurate description of the perturbative Adler function in the energy range Λ 2 QCD Q 2 4 m 2 c . The flavour-diagonal (i = j) correlators contribute to D(Q 2 ) through the so-called nonsinglet topology with the two electromagnetic currents connected by one quark loop. The result can be written as with µ the renormalization scale. Additional disconnected diagrams with each current in a separate quark loop (singlet topology), which are also present for i = j, start to contribute at order α 3 s , but with three massless quarks the flavour trace of both quark loops cancels in the sum (Q u + Q d + Q s = 0), leading to a completely negligible effect of O α 3 s m 4 s Q 4 for D(Q 2 ) once the non-zero masses are taken into account.
The (0) superscript indicates that we have not yet incorporated any quark mass correction. Since the Adler function is independent of the renormalization scale, one can trivially reconstruct the coefficients K n,p with p > 0 from those with p = 0, simply taking into account the Renormalization Group Equation (RGE) satisfied by the strong coupling: The coefficients of the Adler function are known up to five loops, i.e. at O(α 4 s ), while the QCD β function has been already computed to O(α 5 s ). We collect the values of K n≤4,0 , β n≤5 and the rest of perturbative coefficients in App. A. Since we need to truncate the series, a residual scale dependence (of the first unaccounted order) remains. In order to avoid higher-order corrections enhanced by large logarithms of the renormalization scale, one should set µ 2 = ξ 2 Q 2 with ξ 2 a number of order 1. The exact choice is however ambiguous and a priori arbitrary, just as it is the exact scheme choice of how to minimally subtract when renormalizing. Modifying the residual scale dependence through the variation of ξ 2 in a reasonable interval around unity can then be used to estimate perturbative uncertainties [40]. One has In order to numerically evaluate the Adler functions below the charm threshold, we then need α s (µ 2 ) with n f = 3. The standard input is however α s (M 2 Z ) with n f = 5. One can translate one into another by supplementing the RGE given above with the decoupling relations, is the running mass of the heavy quark that has been integrated out and a /π, which is also needed as input. The running masses satisfy the following renormalization group equation where the perturbative γ n coefficients are known up to n = 5. We will take as extra input m c (m 2 c ) = 1.275 (5) GeV and m b (m 2 b ) = 4.171 (20) GeV from the FLAG lattice review [41][42][43][44][45][46][47][48][49][50][51][52][53][54] and perform the decoupling at µ = m q (m 2 q ). Quark-mass uncertainties are small enough to be negligible. Our results for α s (Q 2 ) ≡ α  Table 1 for several choices of α (n f =5) s (M 2 Z ) and Q, at different orders in the strong coupling. We have checked that the α s (Q 2 ) values fully agree with the corresponding values obtained with RUNDEC [55]. In Fig. 2 we show the corresponding Q 2 dependence for α , adding as perturbative uncertainties the quadratic sum of variations due to changing K 5 in a conservative interval (−125, 675) and ξ 2 ∈ (0.5, 2) [36]. The range chosen for K 5 includes the values advocated by renormalon models, Padé approximants, effective charges and conformal mappings [56][57][58][59][60][61][62], but also allows for a correction of opposite sign. The interval of variation for the renormalization scale is conventional for low-energy analyses, e.g. see [63,64]. This is partially justified by the fact that taking a too small value for ξ 2 , one would be using an ill-defined expansion parameter α s (ξ 2 Q 2 ) without any real justification, leading to unreasonable uncertainties. 5 Alternative prescriptions to circumvent this issue, such as asymmetric scale variations, can be found in the literature. See for example Ref. [65]. (Q 2 ) are given in the right columns at different orders in α s (Q 2 ). Our central value for the fifth-order coefficient, K 5 = 275, has been adopted in the last column.

Strange mass corrections
The light quark masses are not exactly zero, and this leaves a small imprint, fully dominated by the strange quark mass. Let us then safely neglect the tiny effects suppressed by O( where m s (Q 2 ) ≡ m (n f =3) s (Q 2 ) and the coefficients are once again shown in App. A. Numerical values for the associated corrections at different Q, taking as input m s (µ 2 0 ) = (92.03 ± 0.88) MeV at µ 0 = 2 GeV [41,42,46,47,49,[66][67][68][69][70], can be found in Table 2. 6 Let us note the very bad behaviour of the perturbative series (3.9), which appears to show its asymptotic behaviour from the first terms. Fortunately, the whole mass correction is very suppressed by the small value of the strange quark mass and its electric charge. We will take as central value the average between truncating at O(α 2 s ) and O(α 3 s ) and half their difference as an additional perturbative uncertainty. The associated Q 2 dependence, together with our estimated error bars are shown in Fig. 3.

Heavy-quark mass corrections
Internal heavy-quark loops induce charm-mass corrections into the light-quark correlators, which are suppressed by powers of Q 2 4m 2 c and start to contribute at O(α 2 s ). In order to take into account these small, but sizable when approaching the charm threshold, contributions to the Euclidean Adler function, we can use the known results at order α 2 s for the associated contributions to R(s). Two distinct topologies appear there. One corresponds to the fourquark cut, ρ R , which starts at s = 4 m 2 c . The second corresponds to the vertex correction ρ V , which, in the chiral limit, starts at s = 0. The exact expressions can be found in [71]. Taking into account that these are the only cuts induced by those topologies, one can  reconstruct the associated contribution to the Adler function by using the same kind of dispersion relation as above: where C F T F = 2/3. As expected, for small enough values of Q 2 both integrals admit expansions in powers of Q 2 4m 2 c . For the former, the Q 2 4m 2 c expansion can only be performed after integration, since s ∈ (0, 4m 2 c ). However, one can first expand ρ V (s) in powers of s 4m 2 c whose leading, next-to-leading and next-to-next-to-leading terms can also be found in Refs. [72,73]: In the right panel of Fig. 4, the agreement between this expansion and the result obtained computing the integral numerically, using the full expression for ρ V (s), can be seen. Following an analogous criteria as for the light quarks, we will estimate uncertainties from higher orders by changing µ between mc(m 2 c ) √ 2 and √ 2 m c (m 2 c ). 7 Similar arguments can be used for the second term. In this case, the expansion in powers of Q 2 can be performed before or after integration. In the first case, the integrals can be computed numerically using the full expressions for ρ V (s) and ρ R (s). On the other hand, to perform the Q 2 expansion after integration, we have expanded ρ V (s) + ρ R (s) first at large s, which agrees with Ref. [71], in order to integrate them analytically. Both ways yield very similar results since, as can be seen in Fig. 5, the expansion of ρ V (s) + ρ R (s) agrees very well with the full expression. The contribution to the Adler function then is: (3.14) 7 This is just one possible choice to circumvent the overestimated uncertainty due to an ill-defined ex- were taken. For heavy quarks an alternative way of achieving this, based on using different scale choices for mc(µ 2 ) and αs(µ 2 ), which is useful to avoid underestimating uncertainties in fits with combined moments, can be found in Refs. [74,75].  Figure 5. Left: Comparison between the exact form of ρ V (s) + ρ R (s) and its expansion in the limit 4m 2 c /s → 0. As can be seen, they are equal for 0 < 4m 2 c /s < 1. Right: Comparison between D L,mc ii,> (Q 2 ) computed numerically and computed analytically, using the expansion (3.13) for ρ V (s)+ρ R (s) and performing an expansion in Q 2 /4m 2 c afterwards. The input value α 1184 has been adopted.
The comparison between this expansion and the exact result is plotted on Fig. 5. As expected, the expansion in powers of Q 2 /4m 2 c diverges from the exact result near threshold. It is worth making a parenthesis to note how this second integral fully dominates if one takes instead Q 2 4m 2 c . In that case one may use the expansion (3.13) of ρ R (s) + ρ V (s) at s m 2 c . The leading term diverges as s → ∞, generating a large α 2 s contribution to D(Q 2 ) that grows logarithmically at large Q 2 : However one should keep in mind that the strong coupling with n f = 3 flavors should not be used as expansion parameter far above the charm threshold. In this limit, at order α 2 s , one then has which shows how the logarithm of the charm mass is properly reabsorbed into the n f = 4 strong coupling, through the QCD matching conditions, while the constant α 2 s term reproduces the known n f dependence of K 2,0 .
The total heavy-quark contribution to the light Adler function, D L,mc , is plotted in Fig. 6. The expanded expression in powers of Q 2 /(4m 2 c ) turns out to provide an excellent approximation to the exact numerical result in the full range of Q 2 values analysed.
The associated Adler function is then (3.20) The C (0,1,2) j coefficients are known up to j = 30 [79,80] and C (3) j up to j = 3 [81][82][83][84]. The contribution to C 4 from topologies associated to quark-connected (non-singlet) contributions has also been computed [85], and very good approximations to the coefficients C (3) j from j = 5 to j = 10 are also known [86,87]. They are typically given at the renormalization scale µ = m c (m 2 c ), but one can then trivially recover them at arbitrary scales by using RGEs. We compile them in App. A. In Fig. 7 we show the associated Adler function up to different orders in α s , cutting the series at j = 10. 8 The series is observed to stabilize after including the two-loop corrections. We also show in Fig. 8 the convergence of the energy expansion at three loops, truncating at several values of j. As expected, the energy series breaks down slightly below Q 2 ∼ 4m 2 c . We find that one actually needs to keep quite high orders. Nevertheless, for practical purposes j = 10 is high enough to safely neglect higher orders below Q 2 ∼ 5.5 GeV 2 . In order to estimate perturbative uncertainties we will vary the renormalization scale in the interval µ = ( 1 . 9 This, together with the uncertainty coming from the input value m c (m 2 c ) = 1.275 (5) GeV, are the main sources of error from this contribution.
Finally, from j = 4 a singlet topology with a massless (three-gluon) cut also appears at four loops, whose leading contribution at Q 2 m 2 c is driven by a known logarithm [88] Π s cc (Q 2 ) = − 17d abc d abc 243000 with d abc d abc = 40/3. While this leads to a logarithmic divergence in the associated coefficient, the limit Q → 0 does not give any problems for the Adler function itself. One finds Based on the known relative values of the constant coefficients in the analogous axialcurrent singlet contributions (starting at α 2 s ) [79], one expects C = 0 ± 3, so the predictive power for this tiny correction is very limited. Taking into account this uncertainty and the one from changing the residual scale dependence of O(α 4 s ) in the interval 8 For the fourth loop and j > 4 we take the approximate coefficients obtained in [87], which succeeded in giving an excellent prediction for the nowadays exactly known j = 4 [85]. 9 The rationale behind this asymmetric choice is that j = 10 is not precise enough near Q 2 ∼ 5.5 GeV 2 if µ = √ 2 mc(m 2 c ) is taken, which has nothing to do with the perturbative uncertainties of the Taylor coefficients that we want to estimate.
1184 has been adopted. . D cc (Q 2 ) at order α 2 s and at different orders in the expansion of Q 2 /4m 2 c . As expected, the convergence is better the higher the order of the expansion is. At the same time, when going beyond the radius of convergence, the higher the order the faster it goes to infinity. loop corrections associated to the bottom quark are incorporated in a completely analogous way. The uncertainty of this contribution is dominated by the error on the input value of the bottom-quark mass.

QED corrections
At the precision level that we have, it is worth assessing the size of the leading QED corrections. For large Euclidean momenta, QED corrections to the Adler function can be computed perturbatively just as one computes the QCD ones.
At this level, there is a subtlety to be considered. Technically, in a full computation of the Adler function in QED plus QCD, one should incorporate the disconnected topology corresponding to a photon between two quark loops. However, this contribution is part of the vacuum polarization of the photon propagator, which, by definition, does not enter into R(s). Then, consistently, at least with R-ratio data, 10 we will not take that contribution into account. The remaining leading QED corrections are well known. Up to O(αα s ) [89] and heavy-quark effects, the relevant contributions can be taken into account by the following shift Taking as central value the average between the results for D L(0) ii obtained at Q = 2 GeV with α s = 0.1184 and α s = 0 (both choices are correct up to O(αα s ) effects) and half their difference as perturbative uncertainty, one finds, where the error includes a conservative estimate of missing QED corrections associated with heavy-quark loops. Perturbative QED corrections are then negligible, at the current precision level.
In Fig. 10 we plot our final result for the perturbative Adler function D(Q 2 ), including all estimated uncertainties. 12 The relative size of the different contributions and their corresponding errors are shown in Fig. 11.

Adler function based on the experimental ratio R(s) and lattice data
In this section we discuss in some detail the experimental data that we use for R(s) (based on the DHMZ compilation [24,28]) and the lattice inputs (based on the Mainz results of 11 The only ones not yet discussed are the 3 i=1 QiD ic(b) ones, which can only enter through disconnected diagrams. However, taking into account that 3 i=1 Qi = 0, they vanish in the chiral limit and then they are suppressed by α 3 , which makes them completely negligible. 12 Let us note how the observed bending at Q 2 ∼ 5.5 − 6 GeV 2 is not a physical feature of the Euclidean Adler function, but a first signature of the breakdown of the series in powers of Q 2mc . As a consequence we will restrict the comparisons with other determinations to Q 2 < 5.5 GeV 2 .  (47)(24)     Ref. [14]).

Experimental ratio R(s)
The experimental results on R(s) used in this work are based on the DHMZ data compilation made in Refs. [24,28]. More concretely, we make use of a previous study of the different contributions to R(s), performed with a full treatment of the uncertainties and their correlations. This information allows us to evaluate the hadronic running of α(Q 2 ), ∆α had (Q 2 ) (see the top panels of Fig. 12), which is then used in order to derive the Adler function with its various uncertainty components and its covariance matrix. The Adler function emerging from the R(s) data is displayed in the bottom panels of Fig. 12, together with the associated correlation matrix. The relative size of the various contributions to D(Q 2 ) and of their corresponding uncertainties are shown in Fig. 13. At low values of Q 2 , the exclusive channels fully dominate. This contribution, corresponding to the region √ s th < √ s < 1.8 GeV, is derived based on the measurements of 32 exclusive channels. This compilation takes into account the statistical and systematic correlations between the different points/bins of a given measurement, between different experiments measuring a given channel, as well as between different channels [24,28,96,97]. In this procedure, possible tensions between different measurements of a given channel are also taken into account. This is generally done in the combination, through an enhancement of the uncertainties by a factor χ 2 /ndof, applied in all the √ s bins where this factor is larger than unity. In addition, an extra uncertainty accounting for the systematic deviations between the BaBar [98,99] and KLOE [100][101][102] measurements in the 2π channel has been included for the first time in Ref. [28], by comparing the combination results obtained when excluding either of the two experiments. 13 This uncertainty turns out to be dominant in 13 Note: Recently, a new precise measurement of the 2π channel performed by the CMD-3 collaboration has been made public [103]. In the dominant ρ-resonance region, it features larger cross-section values compared to all previous experiments, in particular the most precise ones from BaBar [98,99] and KLOE [100][101][102].
the case of the theoretical prediction for the anomalous magnetic moment of the muon, hence the importance of fully taking this systematic effect into account. The remaining data-based contributions come from the 3.7 GeV < √ s < 5 GeV interval, the dispersive integrals being evaluated based on the inclusive measurements available in this range, and from the narrow J/Ψ and Ψ(2S) resonances.
Unfortunately, no precise enough data are yet available for the remaining regions, 1.8 GeV < √ s < 3.7 GeV 14 and √ s > 5 GeV, and there one needs to rely on perturbation theory for R(s). 15 This happens to be a more critical and limiting factor for the associated Adler function, especially at large Q 2 , where this contribution eventually dominates when Q 2 GeV.

Lattice Adler function
The study about the hadronic running of the electromagnetic coupling from lattice QCD presented in Ref. [14] contains all the needed details to extract the corresponding Adler function together with estimated uncertainties and correlations. In order to be able to keep track of them, we use the rational approximation forΠ(Q 2 ) presented in that reference: where x ≡ Q 2 /GeV 2 and the correlation matrix of the expansion coefficients is given by It is worth noticing that, when comparing the rational approximation ofΠ(Q 2 ) with the tables of that reference, we observe that, while this approximation gives a very accurate description of the central value up to Q 2 ∼ 7 GeV 2 , a significant reduction of the uncertainties (up to 50%) starts appearing at 2 GeV 2 , due to a more conservative treatment of discretization effects and, to some extent, to the constraint due to the assumed rational approximation ansatz. We obtain the corresponding lattice Adler function by simply using Eqs. (2.1) and (2.2). The result is displayed in Fig. 14.
While dispersive integrals computed with these new inputs are certainly enhanced and closer to the ones obtained from Lattice QCD, it is of upmost importance to first achieve a better understanding of the tensions on the experimental side. In particular, one needs to understand the source of tension between the CMD-3 results and the former CMD-2 measurements [104,105], performed in somewhat similar conditions by the same group.
14 The relatively precise BES III and KEDR results are in tension in the range 3.40 GeV < √ s < 3.67 GeV [106]. 15 Notice, however, that the perturbative uncertainties include also the estimated size of potential violations of quark-hadron duality in the region 1.8 − 2 GeV [28].

Comparison of the three different approaches to D(Q 2 )
We can finally perform the comparison of the three descriptions of the Adler function.
For the perturbative one we will take the same inputs as above i.e. α  Fig. 15. In Fig. 16 we quantify the tension among the different descriptions of the Adler function in terms of the statistical significance of their differences, i.e.
, (i, j = pQCD, e + e − data, latt). (5.1) We observe the following: 1. The estimate based on e + e − data has smaller uncertainties than the lattice determination. The quoted pQCD precision becomes competitive at 2 − 3 GeV 2 .
2. In the whole energy range analysed, the lattice determination of D(Q 2 ) has a larger central value than the one inferred from e + e − data. This follows the same trend as in g − 2 and ∆α had . However, within the quoted uncertainties, these two estimates of the Adler function remain compatible at large values of Q 2 , the statistical significance of their difference being ∼ 1σ.
3. The lattice determination is in excellent agreement with pQCD in the region where perturbation theory is reliable, i.e. at Q 2 2 − 4 GeV 2 . The differences between the two determinations (central values) is only ∼ 0.2σ. Some tension, larger than 1σ, is observed at lower values of Q 2 . This is expected, both because systematic perturbative uncertainties become less reliable and because power corrections are expected to emerge.

A significant tension between the determinations of D(s)
from e + e − data and pQCD emerges below 5 GeV 2 . It is larger than 2σ, and it even surpasses 3σ at Q 2 ∼ 1.5 − 3 GeV 2 . Let us give three possible explanations for it. First, at low Q 2 values, power corrections can become sizeable and may explain the discrepancy. We analyze this possibility in more detail in the next section. The second possibility is that the strong coupling is lower than the value used as an input. From that perspective one can translate this tension into a tension between the value of α s obtained from a fit to e + e − data and the lattice average that we are using as input. We will study this perspective in more detail below. New precise results from novel methods agreeing with the lattice average [108] suggest that such a large disagreement is unlikely. The third possible explanation for the tension at large Q 2 values, appears to be possible unaccounted systematic effects in the dispersive evaluation of the Adler function (see Fig. 13 for the various contributions to this evaluation and to its uncertainty).

Nonperturbative corrections to the perturbative Adler function
We aim to assess up to which level power corrections can account for the deviations observed in the previous section between the Adler function emerging from e + e − data and pQCD. This is not straightforward because the needed vacuum expectation values are not known from first principles and, in general, their numerical values can depend on the way one truncates the (asymptotic) perturbative series. In fact, the factorial growth of the perturbative expansion at large orders generates infrared ambiguities (when one tries to reconstruct the Adler function from its Borel sum) that scale as inverse powers of Q 2 and are expected to be reabsorbed into the nonperturbative terms of the OPE. One may be tempted to state that those effects are already accounted for in the perturbative systematic uncertainties. However, the existence of vacuum matrix elements is well established beyond perturbation theory [2,3]. The nonperturbative nature of the QCD vacuum generates non-zero vacuum expectation values for many composite operators such as the quark condensate 0|qq|0 , responsible for the breaking of chiral symmetry, or the gluon condensate 0|a s G µν G µν |0 , which breaks the scale invariance of massless QCD. Pure nonperturbative observables, where perturbation theory vanishes, allow for cleaner determinations of the corresponding vacuum condensates because their numerical effects cannot be masked by perturbative uncertainties. This is the case of the two-point function of a left-handed and a right-handed currents (the V V − AA correlator), which is identically zero to all perturbative orders in α s but receives non-zero contributions from D ≥ 6 vacuum condensates that are order parameters of the chiral symmetry breaking. The sizes of the leading power corrections to this correlator are well known, since they can be directly extracted from the τ decay data [109,110]. Since there is no reason to neglect them, neither in the vector correlator nor in the axial one, it is then a must to incorporate power corrections for a complete description of the OPE-based Adler function.

Light-quark correlators
The leading nonperturbative contribution to the Adler function is, up to negligible up and down quark-mass corrections [9,36], The numerical value of the gluon condensate is quite uncertain [111,112], since it is difficult to separate its effect from the ambiguity generated by the asymptotic tail of the perturbative series, which is supposed to be already included in the perturbative uncertainty. On the other hand, from general grounds we know that the gluon condensate is positively defined [2]. This is an important point, because it actually means that its corresponding D = 4 power correction goes into the wrong direction to explain the tension between pQCD and the experimental data on R(s). To be on the conservative side, let us take the central value estimated in Ref. [2], but with a 100% of uncertainty, i.e. α s π GG = (0.012 ± 0.012) GeV 4 . (6. 3) The strange quark condensate is better known because it is related to the kaon mass and decay constant by chiral symmetry [113]. Since it is an order parameter of the chiral symmetry breaking (it vanishes to all orders in perturbation theory), the quark condensate does not suffer from the perturbative ambiguity mentioned before. At lowest-order in chiral perturbation theory it gets determined by the old Gell-Mann-Oakes-Renner relation [114]. However, it receives large higher-order corrections that enhance its final uncertainty [115][116][117]: The combined dimension-four correction to the electromagnetic Adler correlator in Eq. (6.2) takes then the value: Let us also account for the D = 6 contribution. Up to residual pieces that vanish in the electromagnetic sum, one can write the D = 6 contribution as It is convenient to rewrite The O 6,V −A contribution is a genuine vacuum condensate whose nonzero value, O 6,V −A ≈ −0.0035 (9) GeV 6 [109], is well established and understood beyond perturbation theory, and its effect is unrelated to the perturbative series, which is identical for the vector and axial channels. Nonperturbative effects in the observed spectrum (see for example [118,119]) are known to be suppressed for the V + A combination with respect to the V − A, which motivates to assume |O 6,V +A | < |O 6,V −A | [120]. This inequality holds (by far) in the large-N C limit, which gives O ∞ 6,V +A = − 2 9 O ∞ 6,V −A , reproducing the old vacuum saturation approximation, which is also known to work well in predicting O 6,V −A and some rigorous inequalities [2]. Taking this into account, we will adopt Notice that the assigned uncertainties to those corrections potentially contaminated by perturbation theory (gluon condensate and O 6,V +A ) are above 100% of their corresponding estimates, guaranteeing that any potential double-counting effect is consistently absorbed by our conservative errors.

Charm correlator
Much less relevant for our analysis is the contribution of power corrections to the charm correlator. The leading power correction is given by the gluon condensate contribution [2]: Taking again a s GG = (0.012 ± 0.012) GeV 2 , one can easily check that the contribution to D(Q 2 ) remains below 10 −3 and can then be neglected.

Discussion
The nonperturbative correction to the Adler function is then given by a s a s GG + 4 1 + a s 3 + 27 8 + 4 ζ 3 a 2 s m s ss ] + 24π 2 O 6,V Q 6 . (6.12) Adopting the conservative numerical estimates in Eqs. (6.3), (6.4) and (6.8), we find that the nonperturbative uncertainty is actually larger than the perturbative one at low Q 2 values and, at the current precision level, it cannot be fully neglected even at Q 2 ≈ 4 GeV 2 . We show in Fig. 17 the comparison of the OPE (pQCD plus condensates) Adler function with the other approaches, and in Fig. 18 the statistical significance of the differences with respect to the lattice and e + e − -data evaluations of the Adler function. We observe that the origin of the tension at ∼ 1 GeV 2 between the perturbative prediction and both lattice and DHMZ results can indeed be explained by genuine nonperturbative effects. Incidentally, the input central values assumed for the vacuum condensates fit very well the shape of the distribution, although slightly different estimates cannot be discarded within the current experimental and lattice uncertainties, also depending on the size of higher-order power corrections. On the other hand, for Q 2 2 GeV 2 , the observed tension between the analytic Adler function and the data-based one gets slightly reduced. This reduction decreases with Q 2 , as the effect of the nonperturbative terms diminish, and a tension of up to ∼ 2 σ remains.

Determination of α s from the Adler function
Instead of using α s as input to compare the perturbative D(Q 2 ) with the other approaches, we can reconvert the comparison into an α s extraction. The extraction can be done at each value of Q 2 by solving the equation: for α s while keeping Q 2 fixed. D OPE (Q 2 , α s ) is the sum of the perturbative and nonperturbative contributions of the theoretical Adler function, Eqs. (3.26) and (6.12), keeping α s as a variable. D data (Q 2 ) is the Adler function obtained through either the experimental ratio R(s), D R(s) (Q 2 ), or the lattice data, D latt (Q 2 ), which were discussed in Sec. 4. In this section we illustrate the procedure by using the Adler function based on R(s). The results for the analogous lattice fits are relegated to App. C. 16 16 Notice however how lattice data is not limited to the full EM correlator and then better strategies to extract αs can in principle be pursued. As explained in Sec. 4, the R(s) data relies on perturbation theory for the regions 1.8 GeV < √ s < 3.7 GeV and √ s > 5 GeV. This means that the R(s)-based Adler function also contains a residual dependence on α s , which we take into account in the extraction. We decompose D R(s) into the sum of two contributions: one coming from experimentally measured values, D R(s) exp (Q 2 ) (see Fig. 13), and one corresponding with the perturbation theory contribution, D R(s) P (Q 2 , α s ). The resulting equation to solve for α s at each Q 2 becomes: The values for α (n f =5) s (M 2 Z ) obtained from both lattice and R(s) data, together with their uncertainties, including both all experimental and theoretical sources, are plotted in Fig. 19, as a function of the Q 2 value at which they are derived. As expected, the same observations made when comparing the different Adler functions are applicable here as well. That is, 1. The central values of the strong coupling extracted from e + e − data are smaller than the ones extracted using the lattice determination of D(Q 2 ).
2. The values of α s extracted from lattice data are in agreement with the FLAG lattice average.
3. The α s extracted from e + e − data is between 1.5 and 2 σ below the FLAG lattice average in the 3 GeV 2 < Q 2 < 5 GeV 2 region. Notice however how the different Q 2 points are strongly correlated among each other.
In the following, we discuss how these different values for α s can be compared quantitatively and combined.

Averages
A potential improvement in terms of precision for the extracted α (n f =5) s (M 2 Z ) could come from combining the values obtained at different Q 2 values, taking into account their correlations. One possibility would be to perform the combination through the minimisation of a χ 2 function, with respect to α av s . In this function,ᾱ extr s is the vector containing the extracted values of α s at each Q 2 , C is their covariance matrix 17 andᾱ av s is a vector with the same dimension asᾱ extr s containing the parameter α av s . However, one is faced with the following limitations: 1. As discussed in the previous section, the nonperturbative contributions to the Adler function are a large source of systematic uncertainty at low Q 2 , which additionally is not fully controlled. This can clearly be seen in Fig. 19. Therefore, one should avoid using the extracted values in this energy region.
2. The values extracted from e + e − data display strong experimental correlations (see Fig. 12). These correlations are stronger between neighboring points and a larger precision in the estimation of the corresponding covariance matrix would be needed in order to find meaningful results for the combination, which may otherwise not be realistic. Indeed, if one includes too many consecutive points, eventually one is going to find a covariance matrix with null eigenvalues, which would imply absolute predictions (i.e. no uncertainties in certain linear combinations). Actually, this is a consequence of the approximations involved in finding the original covariance matrix. 18 Additionally this makes the experimental correlation matrix singular, which therefore cannot be inverted. As a result, a χ 2 function cannot be constructed for the whole set of extracted values.

3.
A similar problem arises for the theoretical correlation matrix. In general one expects that the perturbative uncertainties are dominated by the first unknown coefficient, in this case K 5 , and the nonperturbative ones by the first unknown power correction, say O 6 . We have supplemented the perturbative uncertainty by renormalizationscale variations, which in general one expects to account for uncertainties due to higher-order effects. However this is going to fail if one artificially looks for linear combinations of data points such that either K 5 , O 6 and/or the scale variations cancel or they appear suppressed, for example by numerical prefactors, with respect 17 In principle, for e.g. relative uncertainties, this covariance matrix can itself depend on α av s , which could be addressed through an iterative fitting approach (see e.g. Refs. [121,122] and references therein) or using fitted nuisance parameters applied as (constrained) scaling factors of the theoretical prediction [123,124]. However, here the input α extr s values are similar (see e.g. Fig. 19) and the impact of this effect is small. 18 An illustrative example of this (extreme) case, this time for fit inputs from lattice QCD, would consist in simply starting from Eq. (4.1) and taking 7 different points. While the rational approximation and their associated uncertainties are expected to work generally well both inΠ(Q 2 ) and linear combinations of it, it would clearly fail in predicting the uncertainties of the linear combination corresponding to the zero eigenvalue (corresponding to a null estimated uncertainty), which would directly dominate (because of an apparently infinite precision) the determination of any theoretical parameter depending on it.
to the contributions of higher-order coefficients, which then by construction are not going to be negligible with respect to the accounted effects. In a naive χ 2 fit, the extracted value of α s is going to be dominated by the most precise linear combination of data points, taking into account those theoretical uncertainties, which are precisely the directions where the estimators are prone to underestimate them, leading to very aggressive predictions. An explicit example of this kind of direction is the logarithmic derivative of the Adler function, implicit in fits to consecutive data points. 19 Indeed by taking consecutive (adimensionalized) derivatives one is going to trigger the sooner breakdown of the OPE. Schematically, for the nonperturbative contributions to D, D NP , one has 20 and then higher-dimensional corrections become more and more important with respect to lower-dimensional ones at a fixed energy, eventually leading to potentially underestimated theoretical uncertainties in those directions. Analogously, the scale variation is going to be in general a good estimator of perturbative uncertainties, because its variation is in general of the same order as the neglected perturbative contributions. However, one would clearly fall into one version of the well-known look-elsewhere effect if one artificially looks for directions in which it happens to cancel: the fact that the scale-dependence accidentally cancels in some linear combination does not guarantee that the contributions from higher-order coefficients are also cancelling. There is also the possibility that new topologies emerge at higher orders, inducing an uncertainty not well accounted by the scale variations.
In summary, the χ 2 function defined in Eq. (7.3) and the result of its minimisation are sensitive to uncertainties on the uncertainties and on the correlations (i.e. to uncertainties on the covariance matrix), present for both the experimental and the theoretical components. Starting from remarks made in the context of ATLAS jet performance and cross-section studies [125][126][127], the relevance of the uncertainties on the covariance matrices (in particular for what concerns the implications for combination methods) has been pointed out in the context of the theoretical predictions for the anomalous magnetic of the muon [28,128,129]. More recently, similar remarks about the uncertainties on uncertainties were made for what concerns the procedure of quantifying the significance of the data-theory tensions, in this same context [130].
Taking all these aspects into account, we will restrict to data sets with at most three points and Q 2 3 GeV 2 . We will consider the two sets of three Q 2 points shown in Table 4. The first set, set 1, has a better behaviour in the expansion in powers of Q 2 /m 2 c , whereas 19 In fact we find that the eigenvectors associated to the lowest eigenvalues of the experimental covariance matrix correspond, in first approximation, to the highest-order derivatives in the discrete approximation. They can be related to different, more localized, weights when integrating R(s). 20 An analogous issue occurs for the Q 2    Table 4. The two sets of three Q 2 points chosen for the averages with their experimental (exp), perturbative (pert) and theoretical (th) symmetrized uncertainties. By perturbative uncertainty we mean the uncertainty coming from the use of perturbation theory in the regions 1.8 GeV < √ s < 3.7 GeV and √ s > 5 GeV for the ratio R(s).
the second set, set 2, is less affected by potential nonperturbative effects. Additionally we will also consider two variations with only two Q 2 values, set 1 * and set 2 * , where the midpoint of the corresponding set has been removed. The only remaining input necessary to compute the χ 2 is the covariance matrix C. In order to evaluate it, we can use linear error propagation, for both the experimental and theoretical components of this matrix. In particular, for the theoretical component this is justified, since the Adler function is an approximately linear function of α s in the energy region we are considering. However, there are different possible choices aimed to avoid the issues mentioned above when assigning correlations between the theoretical uncertainties at different points. We will explore the different possibilities in Sec. 7.2.
When performing such combination, we must also take into account the possibility of the χ 2 minimisation yielding a biased result, caused by the limited precision with which the covariance matrix is known. According to the Gauss-Markov Theorem, the minimisation of the χ 2 in Eq. (7.3) is equivalent to performing a weighted average, while optimizing the weights (constraint to have the sum equal unity) such that the uncertainty of the average is minimum [131,132]. This yields with1 being a vector of the same dimension asᾱ s extr , with all its entries equal to 1. As a consequence, if the covariance matrix is poorly known, these weights can be biased (in some cases they can e.g. take negative values or values larger than unity, which could make the average value to be outside the range of the extracted values). Therefore, we have to consider alternative averaging procedures that, although will yield results with somewhat larger uncertainties (see Sec. 7.3), are free of any bias caused by implicit assumptions in the derivation of their weights. We have considered the following: 1. A simple average, where all the inputs have the same weight. This is, the inverse of the number of input values considered.

2.
A weighted average where the weights are proportional to the inverse of the experimental uncertainty squared. 21 The averaged values obtained using these different averaging procedures for the four sets considered can be found in Sec. 7.3.

Theoretical uncertainties
As discussed throughout the text, there are different sources of theoretical uncertainty to the Adler function and for some of them assigning a 100% correlation between different Q 2 values may not be the best choice. This is for example the case of the scale variations (it is probably not realistic to assume that there is a single renormalization-scale choice such that the perturbative uncertainties of all truncated series can be removed) or the nonperturbative one parameterized by O 6 (knowledge of O 6 would not completely remove the nonperturbative uncertainty, since we would need to account for the higher-dimensional contributions). We will then repeat the same fits under four assumptions on the correlations between the theoretical uncertainties at different Q 2 : 1. All correlations for the same "theoretical source", renormalization scales or O 6 , are 100%.
2. All correlations for the same "theoretical source" are 100% except for the scale variations, that are assumed to be uncorrelated among points with large separations in Q 2 .
3. All correlations for the same "theoretical source" are 100% except for the O 6 one, which is assumed to be uncorrelated among points with large separations in Q 2 .
4. All correlations for the same "theoretical source" are 100% except for the scale uncertainties and the O 6 ones, which are assumed to be uncorrelated among points with large separations in Q 2 .

Results
The numerical results obtained with the four different assumptions on the theoretical correlations are displayed in Tables 5, 6, 7  (M 2 Z ) extracted from the four different choices of Q 2 points (sets 1, 2, 1* and 2*) and with the three different averaging procedures: χ 2 minimisation, simple average and weighted average. Clearly, the results in Table 5 have too large χ 2 values associated to them, both for sets 1 and 2. As explained above, this is not necessarily a signature of inconsistent experimental data sets, but can also be due to the limitations explained in Sec. 7.1. Much more reasonable values are obtained by taking points with broader separation in Q 2 , i.e. the set 1 * , or with the alternative choices for the theoretical correlations proposed in Sec. 7.2. 21 In Ref. [133] one can find further discussions on unbiased combination procedures. These feature realistic uncertainty estimates, based on the propagation of the full information on the uncertainties of the inputs, with their correlations.   All in all we find that there is not much gain in combining several Q 2 values. We observe that the minimisation procedure yields a bias towards larger values of α s , which are preferred by the (potentially dangerous) directions associated to small eigenvalues of the original covariance matrix. Once these are treated more conservatively (i.e. Tables  6, 7  exhibiting again the discrepancy between the dispersive and the lattice-based results. Nevertheless, this also shows that, once the situation with respect to the different tensions related to R(s) is clarified, a determination of α    at some reference Q 2 point (in this case 5.41 GeV 2 , yielding the most precise α (n f =5) s (M 2 Z ) value among the Q 2 points considered here) and each of the other considered Q 2 points respectively. We also compute the corresponding relative differences normalised with respect to the former (reference) α Fig. 20). The uncertainties on these differences are evaluated through a linear error propagation, taking into account the full information on the correlations among the uncertainties of the two corresponding α (n f =5) s (M 2 Z ) input values. The significance of the deviation from zero for each of these differences, computed as the difference divided by its uncertainty, represents a test of the RGE within the Q 2 range of the two points that are being considered. 23 In addition, the uncertainty of each such difference provides a measure of the precision within which the RGE test has been performed. This represents an important information, in addition to the Q 2 range and to the outcome of the test itself. Indeed, we observe that the RGE test here is performed within a precision between 1.4 permil and about two percent, depending on the Q 2 range and the correlation assumptions employed for the theory uncertainties. This can be seen in the plots at the bottom of Fig. 20, where we display the differences between the extracted α  Table 4, normalised by the former. 24 23 The significance squared for the simple differences corresponds to the minimum of the χ 2 in Eq. (7.3) (see e.g. Ref. [134]). Very similar values are obtained for the significance of the relative differences squared (up to small non-linear effects in the uncertainty propagation). We have also checked this correspondence numerically in the current study. 24 The analogous plots using instead lattice data are once again relegated to App. C.  Table 4, together with their uncertainties computed taking into account the correlations. The bottom plots show the same differences normalised with respect to the α The red circle indicates the reference point at Q 2 = 5.41 GeV 2 . For the left plots we assume the theory uncertainties originating from the same source to be fully correlated, whereas for the right plots we assume the theory uncertainties from scale variations and O 6 to be uncorrelated.

Conclusions
In this work we have carefully analysed the QCD predictions for the electromagnetic Euclidean Adler function, below the charm threshold. It was already known that this observable could be studied within pQCD from relatively low energies. From the phenomenological point of view, this is relevant because one can independently determine it at relatively low energies from recent precise e + e − and lattice data, which are known to produce results in clear tension with each other. A comparison between these three different determinations of the Adler function is then compelling because 1. pQCD may discriminate, up to a certain extent, between the e + e − -based and the latticed-based Adler functions if they do not agree. 2. Assuming the validity of one of them, one can test the validity domain of pQCD, just at the edge of the perturbative breakdown. 3. One can study the sensitivity to the QCD coupling of a direct comparison between pQCD and the ratio R(s).
In order to unlock the full potential from known pQCD results at relatively low energies, a consistent analysis beyond the available O(α 2 s ) precision in mass-dependent (decoupling) schemes such as MOM was preferred. Most of this work has consisted in assembling all the needed pieces in an EFT-based MS set-up, to perform such a comparison at O(α 4 s ), with a careful treatment of all associated expansions and uncertainties. This has been done in Sec. 3, where many details are given in order to simplify the reproducibility of our results and facilitate future applications and improvements. The extraction of alternative Adler functions, based on R(s) data and lattice results, has been studied in Sec. 4. We have first compared in Sec. 5 the pure perturbative predictions with the results obtained from both e + e − data and lattice QCD, confirming (although at a somewhat reduced level) the tension observed in g − 2 between these two different approaches. Taking as input the value of α s from the FLAG compilation, α (n f =5) s (M 2 Z ) = 0.1184 ± 0.0008, the pQCD prediction of the Adler function turns out to be in excellent agreement with the lattice determination at Q 2 2 GeV 2 , while the determination from e + e − data lies systematically below it in the energy range Q 2 ∈ [1, 5.5] GeV 2 .
In Sec. 6 we have incorporated the leading nonperturbative power corrections, in order to have a better assessment of the theoretical uncertainties in the lowest range of Q 2 values. Once these corrections are taken into account, the agreement between the QCD OPE prediction and the lattice result extends to the whole analysed range of Q 2 , although the theoretical uncertainties turn out to be large below 2 GeV 2 .
Assuming the validity of the OPE at Q ≈ 2 GeV, the Adler function extracted from the e + e − data lies approximately 2 σ below the OPE predictions. This appears to follow the same trend as the deficit observed in the muon g − 2 integral, when comparing the dispersive e + e − estimate to both the experimental measurement of the muon anomalous magnetic moment and the lattice results, as well as the deficit observed by different lattice groups for the window integral and the deficit in the hadronic running of α when compared to the lattice result from Ref. [14].
In order to fit the Adler function extracted from e + e − data, one would need values of the strong coupling significantly below the current lattice (and PDG) average and/or much larger nonperturbative corrections that would not scale with the expected power corrections. Both possibilities are disfavoured by the lattice data that match beautifully the pQCD predictions, at the achieved precision, even at energies as low as Q ∼ 1.25 GeV.
On the other hand, several possible applications can be expected from our study. Combining the results of this work with future lattice studies could help in assessing the validity domain of pQCD for the different involved correlators at higher precision, while at large Q 2 values, where pQCD is more reliable, the corresponding results can help in understanding discretization effects in the lattice. We also leave for future work exploring an alternative implementation of the Euclidean running of α, based on EFTs (using the MS scheme) supplemented by appropriate energy expansions and possibly interpolations, just starting from the massless descriptions at the different number of flavors, as we have done in this first work at n f = 3. This can also be combined with lattice data, as shown in [14,39].
Once the current tensions between the lattice QCD and data-driven approaches (which are also very relevant in the context of the anomalous magnetic moment of the muon) are resolved, the comparison between data and QCD could provide a determination of α (n f =5) s (M 2 Z ) at the per-cent level. With this respect, forthcoming precise measurements of the hadronic production cross-sections, as well as independent lattice calculations of similar precision, will play a major role. (A.12) Following the notation of Ref. [36], the strange mass corrections are given by a linear combination of three coefficients In general they are all known up to three loops [160][161][162][163],

A.3 Heavy-quark loop coefficients
The contribution to the Adler function from heavy-quark loops depends on the C (n) j (µ) coefficients defined in Eqs. (3.17) and (3.18). In Table 9 we compile the known (nonsinglet) coefficients C (n) j (m c (m c )), up to j = 10 and n = 3, taken from the references referred in the main text. Let us remark that for n = 3, j > 4 they are only approximated values. The n = 3 singlet contribution is given separately in Eq. (3.21).
Using the scale-invariance of the Adler function and the RGE of the running coupling and masses, it is straightforward to find the values of these coefficients at any other renormalization scale. For example, for the four-loop coefficients one finds: 3 (m c (m c ))+ 16544 8505 L 3 + 301549372 13395375 L 2 + 588425644445059 240045120000  Table 9. Non-singlet heavy-quark coefficients C  B Interplay ofΠ 08 with perturbative QCD In Ref. [14], the rational approximation of another correlator, obtained using lattice data, is presented, Π 08 , which is defined as, From Eq. (B.2) it is straightforward to see that the correlator vanishes at all orders in massless perturbative QCD, since in that limitΠ ll =Π 33 and the disconnected topology associated to the second term, which in general starts at O(α 3 s ), vanishes for the sum. The leading OPE contribution toΠ 08 (Q 2 ) is then given by the perturbative mass correction, scaling as ∼ m 2 s Q 2 , so thatΠ 08 (Q 2 → ∞) = Π 08 (Q 2 → 0). In fact, let us note that this asymptotic scaling implies that the correlator satisfies an unsubtracted dispersion relation, The associated (massive) perturbative Adler function is, safely neglecting m u,d , simply given by where ∆ ms D L 33 (Q 2 ) and ∆ ms D 3l (Q 2 ) refer to the strange mass correction of the associated correlators, being ∆ ms D 3l (Q 2 ) = Since the perturbative QCD contributions are suppressed by two powers of the energy, one may expect nonperturbative effects to be numerically more relevant. They only enter suppressed by two extra powers of the energy and, additionally, they are linear instead of quadratic in the small strange quark mass, since the chirality-conserving nature of the vector current insertions can be recovered by combining a chirality-flipping insertion of the strange quark mass with a second one from the quark condensate. The associated contribution is Let us then take the corresponding rational approximation given in Ref. [14], Π 08 (Q 2 ) = 0.0217 (11) x + 0.0151 (12)  and compare it to the perturbative and OPE descriptions. The result is displayed in Fig. 21. In spite of the extremely bad behaviour of the associated pQCD series, which fortunately plays a very marginal quantitative role in the EM Adler function, an excellent agreement to pQCD (with a somewhat arbitrary truncation criteria) appears to emerge up to Q 2 ∼ 5 GeV 2 . However, when incorporating the D = 4 power correction, a slight tension emerges. Given the different scaling both in energy and in strange quark mass of the pQCD series, ∼ m 2 s Q 2 , with respect to the leading power corrections ∼ C Fit results for α s using lattice data For completeness we provide in this appendix the corresponding lattice data results for Tables 4-8 in Tables 10-14. The lattice version of Fig. 20 is shown in Fig. 22. 26 Once again, the very poor behaviour of this perturbative series at small Q 2 can be seen in Table 2.