Novel Method to Reliably Determine the QCD Coupling from $R_{\rm uds}$ Measurements and its effects to Muon $g-2$ and $\alpha(M_Z^2)$ within the Tau-Charm Energy Region

We present a novel method for precisely determining the QCD running coupling from $R_{\rm uds}$ measurements in electron-positron annihilation. When calculating the fixed-order perturbative QCD (pQCD) approximant of $R_{\rm uds}$, its effective coupling constant $\alpha_s(Q_*^2)$ is determined by using the principle of maximum conformality, a systematic scale-setting method for gauge theories, whose resultant pQCD series satisfies all the requirements of renormalization group. Contribution due to the uncalculated higher-order (UHO) terms is estimated by using the Bayesian analysis. Using $R_{\rm uds}$ data measured by the KEDR detector at $22$ centre-of-mass energies between $1.84$ GeV and $3.72$ GeV, we obtain $\alpha_s(M_Z^2)=0.1227^{+0.0117}_{-0.0132}({\rm exp.})\pm0.0016({\rm the.})$, where the theoretical uncertainty (the.) is negligible compared to the experimental one (exp.). Numerical analyses confirm that the new method for calculating $R_{\rm uds}$ removes conventional renormalization scale ambiguity, and the residual scale dependence due to the UHO-terms will also be highly suppressed due to a more convergent pQCD series. This leads to a significant stabilization of the perturbative series, and a significant reduction of theoretical uncertainty. It thus provides a reliable theoretical basis for precise determination of the QCD running coupling from $R_{\rm uds}$ measurements at future Tau-Charm Facility. It can also be applied for the precise determination of the hadronic contributions to muon $g-2$ and QED coupling $\alpha(M_Z^2)$ within the tau-charm energy range.


Introduction
Quantum chromodynamics (QCD) is the fundamental non-Abelian gauge theory of strong interactions.Its running coupling (α s ) sets the strength of strong interaction among quarks and gluons, which is crucial and deserves the best possible precision.The strong running coupling becomes weak at short distances due to the property of asymptotic freedom [1,2], allowing perturbative calculation of physical observables involving large momentum transfer.The strong running coupling in itself is not a physical observable, but rather a quantity defined in the context of perturbation theory, which enters into perturbative QCD (pQCD) predictions for experimentally measurable observables.Its value must be inferred from such measurements and is subject to experimental and theoretical uncertainties [3][4][5][6].
Total hadronic e + e − annihilation rate R is a fundamental observable in QCD, which provides one of the cleanest platforms for determining α s [7].The R value also contributes to the standard model (SM) prediction for the muon anomalous magnetic moment a µ = (g − 2) µ /2 and the QED running coupling evaluated at the Z pole, α(M 2 Z ), e.g., see Refs.[8,9].Till now, many experimental groups have measured the R value.The recent data [10,11] were given by the BES III detector at BEPC II [12] and the KEDR detector at the VEPP-4M e + e − collider [13].A collection of all available R data is given in Ref. [3].Theoretically, the R value has been evaluated in massless pQCD [7,14], and its QCD corrections have now been calculated in the MS-scheme to order α 4 s [15][16][17][18][19][20][21][22].It has been found that the power suppressed finite-quark-mass effects are well under control [23][24][25][26][27][28][29] and the same applies to mixed QCD and electroweak corrections [30,31].
The R for the continuum light hadron (containing u, d and s quarks) production, denoted by R uds , is usually adopted to test the validity of pQCD calculation in relatively low energy region [32,33].The measured R uds (s) excludes the contribution from resonances and reflects the lowest order cross section for the inclusive light hadronic event production through one photon annihilation of e + e − .So it can be directly compared with the pQCD prediction and be directly used to extract α s , or equivalently the QCD scale parameter Λ.
At present the pQCD calculations for R uds are usually analyzed by using conventional scale-setting method, i.e., one calculates the central value by simply setting the renormalization scale µ r equal to the centre-of-mass energy µ r = √ s; and the theoretical uncertainties are estimated by varying the renormalization scale over an arbitrary range such as √ s/2 < µ r < 2 √ s.This leads to conventional renormalization scheme and scale ambiguities, and makes the scale uncertainty one of the most important systematic errors for pQCD predictions.
The principle of maximum conformality (PMC) [34][35][36][37][38] has been proposed to eliminate the conventional renormalization scale-and-scheme ambiguities.The conventional scaleand-scheme ambiguities are caused by the mismatching of the strong coupling and its corresponding coefficients, since its scale is set by guessing.It is noted that the α s -running behavior is governed by the renormalization group equation (RGE), and then the {β i }terms emerged in perturbative series can be inversely adopted for fixing the correct value of α s .The PMC single-scale-setting approach (PMCs) [39,40] determines an overall effective α s (its argument is called as the PMC scale) for any fixed order prediction with the help of RGE.The PMC scale can be treated as the effective momentum flow of the process.It has been shown that the PMC prediction is free of conventional scale ambiguity [41,42], being consistent with the fundamental renormalization group approaches [43][44][45][46] and the self-consistency requirements of the renormalization group [47,48].However there is still residual scale dependence due to uncalculated higher-order (UHO) terms, which will be highly suppressed due to more convergent pQCD series [49].The PMC reduces in the Abelian limit to the Gell-Mann-Low method [50] and it provides a systematic way to extend the well-known Brodsky-Lepage-Mackenzie (BLM) method [51] to all orders.
In this work, we will first adopt the PMCs approach to deal with the perturbative series of R uds (s).Contributions due to the uncalculated higher-order terms will be estimated by using the Bayesian analysis.Then by using the predicted R uds (s) as the basic input, we will extract the value of α s from the KEDR data on R uds (s) and calculate its effect to Muon g − 2 and α(M 2 Z ).The remaining parts of the paper are organized as follows.In Sec. 2, we will present useful formulas for R uds (s) and give a brief description on the calculation procedures.In Sec. 3, we will show the numerical results and discussions.Sec. 4 is reserved as a summary.

Calculation technology
Total hadronic e + e − annihilation rate R(s) is related to the theoretically calculable Adler function D as follows [52], Here the Adler function D is defined as the logarithmic derivative of the hadronic vacuum polarization function Π, which can be written in terms of Π and the photon field anomalous dimension, γ, e.g.[21] γ and Π are given by the perturbative expansions, where d R = N c is the dimension of the quark representation of the colour gauge group.The coefficients , where the superscripts "ns" and "si" denote the non-singlet and the singlet components, respectively.The singlet contribution starts from order-α 3 s , i.e., Π si 0 = Π si 1 = Π si 2 = 0, γ si 0 = γ si 1 = γ si 2 = 0.All these perturbative coefficients γ ns i , γ si i , Π ns i and Π si i up to four-loop QCD corrections can be found in Ref. [21].
Using the perturbative expansions of γ and Π, one then obtains the perturbative expansion for R(s).As for R uds (s), its perturbative expression reads where specifies the known loop level of the QCD correction, the renormalization scale is set to µ 2 r = s.The results for generic values of µ r can be easily recovered by using the standard RGE evolution.The perturbative coefficients r i can be divided into conformal parts (r i,0 ) and non-conformal parts (proportional to β i ), i.e. r i = r i,0 + O({β i }).The {β i }-pattern at different orders exhibits special degeneracies [36,38,53], which lead to ) )

• • •
where It is noted that for R uds , only u, d and s quarks are produced, thus the number of active flavours is n f = 3.Since f =u,d,s q f = 0, the singlet contribution vanishes in the present considered three-flavor case.The anomalous dimension γ also contains n f terms, but it governs the QCD-induced corrections to the running of inverse QED coupling constant α −1 [21], i.e., d(α −1 )/d ln µ 2 r = −α −2 dα/d ln µ 2 r = −4πγ(α s ), and is independent to the running of QCD coupling constant, thus its coefficients γ ns i are kept as conformal coefficients that represent the intrinsic perturbative nature of R(s).Starting from r 3 , terms proportional to π 2 arise due to continuation of the spacelike perturbative results into the timelike domain.These "π 2 -terms" are also called "kinematical terms", and can be predicted from those of lower order.It is necessary to emphasize that, Eq.(2.3)only partially retains the effects due to continuation of the spacelike perturbative results into the timelike domain, and has certain shortcomings (see, e.g., [54][55][56][57][58][59]). As shown in Eq.(2.8), all "π 2 -terms" are nonconformal, thus will be resummed to a certain level in the PMCs scale-setting procedure.
Following the standard procedure of the PMCs approach [39,42], the overall renormalization scale can be determined by requiring all the nonconformal {β i }-terms vanish, the pQCD approximant (2.3) then changes to the following conformal series, where the PMC scale Q * is of perturbative nature and can be fixed up to N ( −2) LL-accuracy, i.e. ln Q 2 * /s can be expanded as a power series over where the coefficients S i (i = 0, 1, 2) read, Eq.(2.10) shows that the logarithmic form ln Q 2 * /s is a power series in α s , which resums all the known {β i }-terms via the RGE, and is independent of µ r at any fixed order.
The resulting conformal series (2.9) with an overall α s (Q * ) provides not only precise prediction for the known fixed-order pQCD series, but also a reliable basis for estimating the contributions from the uncalculated higher-order (UHO) terms.As an estimation of the UHO terms of the perturbative series, we adopt a Bayesian-based approach (BA) [60,61] to quantify it in terms of a probability distribution.The conditional probability density function (p.d.f.) where c (> 0) is a common boundary for the absolute values of all the known coefficients {c l , . . ., c k } and the unknown coefficient c n one wants to evaluate.h 0 (c n |c) is the conditional p.d.f. of c n given c.The conditional p.d.f. of c given coefficients , can be determined by applying the Bayes' theorem, where h(c l , • • • , c k |c) is the likelihood function for c ; i.e., the joint p.d.f. for the coefficients viewed as a function of c, evaluated with coefficients actually obtained in the calculation.The function g 0 (c) is the prior p.d.f. for c.Both g 0 (c) and h 0 (c i |c) depend on the model assumption.Here we use the CH model [60], which suggests: both ln c and c i are equally probable for all their possible values; all the coefficients that we know and that we want to evaluate are mutually independent with the exception for the common bound (c), which results in Using the CH model, we obtain a symmetric posterior distribution for negative and positive c n : a central plateau with suppressed tails [60,61].The knowledge of p.d.f.f c (c n |c l , c l+1 , . . ., c k ) allows one to calculate the degreeof-belief (DoB) that the value of c n belongs to some credible interval (CI).The symmetric smallest CI of fixed p% DoB for c n is, where the boundary c n is defined implicitly by The expression of c n can be found in Ref. [61].We take p% = 95.5% in the following calculation.

Numerical results and discussions
In numerical calculation, to be consistent we shall adopt the -loop α s -running to obtain numerical predictions for R ( ) uds (s).

Basic properties of the pQCD approximant for R uds (s).
The PMC prediction for R uds (s) up to order α 4 s reads, where Both the PMC conformal series (3.1) and the PMC scale (3.2) are scale-independent, which will have residual scale dependence due to uncalculated terms [49].
As a comparison, we present the conventional prediction for R uds (s) by taking µ r = √ s, The UHO coefficients predicted by using BA are r 5,0 /π 5 ∈ [−0.4622, 0.4622] for the pQCD approximant (3.1) and S 3 /π 3 ∈ [−4.4159, 4.4159] for the PMC scale (3.2).More coefficients at lower order predicted by using BA are given in Tables 1 and 2, where the conventional coefficients r i with fixed µ r = √ s are also presented as a comparison.It is noted that almost all the exact values of these coefficients lie within the predicted 95.5% credible intervals.There are only one exception for the conventional coefficient r 4 (µ r = √ s), i.e., the exact value of r 4 (µ r = √ s) is outside the region of the 95.5% credible interval predicted by using the BA based on the known coefficients r i (µ r = √ s) (i = 1, 2, 3).Because the known coefficients of the conventional pQCD series (3.3) are scale-dependent at every order, the BA can only be applied after one specifies the choices for the renormalization scale, thus introducing extra uncertainties for the BA.Such extra uncertainty can be simply evaluated by varying the renormalization scale µ r in some range, such as, √ s/2 < µ r < 2 √ s.This variation range will be labelled as ∆µ r in the following.There are improved models based on the Bayesian analysis, i.e., the geometric model [62] and the abc model [63], which can be applied to deal with the conventional pQCD series with conventional scale dependence ∆µ r , or the PMC series with residual scale dependence ∆Q * .We stress that in Refs.[62,63] ways to deal with the scale dependence within the Bayesian approach have been introduced, and further investigation is left to future work.
MS , represents the Λ parameter of MS scheme in three-flavor QCD.In Table 3, the central values are calculated according to Eq.(3.2) truncated at corresponding accuracy, and the errors ∆Q * are determined by taking the UHO coefficients of ln Q 2 * /s presented in Table 1.To show the residual scale dependence of the PMC predictions, we present R    4) is presented in Figure 2. Figure 1 shows a reduction of the residual scale dependence for the PMC predictions when increasing the order.While the conventional scale dependence of the conventional predictions is moderate when more-and-more loop corrections have been added as show by Figure 2.   uds (s) has been markedly improved after applying the PMCs scale-setting procedure.The conventional prediction (3.3), whose validity range has been demonstrated to be strictly limited to √ s/Λ > exp(π/2) 4.81 [57][58][59], converges rather slowly when the centre-of-mass energy √ s approaches this value, and the corresponding curves start to swerve quite above the boundary of its convergence range.A partial enlarged description for the PMCs predictions R ( ) uds (s)| PMCs with -loop ( = 2, 3, 4) QCD corrections is presented in Figure 4, where the error band of each curve represents the uncertainty from the residual scale dependence.For definiteness, we present the numerical results of R ( ) uds (s) by taking various QCD corrections ( = 2, 3, 4) with three different input √ s/Λ = 5, 8, 12 in Table 4, where the first errors are from the conventional scale dependence ∆µ r of conventional predictions (Conv.), or the residual scale dependence ∆Q * of PMCs predictions (PMCs), and the second errors are from the UHO of the pQCD approximant R ( ) uds (s).Table 4 shows that, the theoretical uncertainty of R ( ) uds (s) is dominated by the scale error, i.e., the conventional ∆µ r for conventional predictions, or the residual ∆Q * for PMCs predictions.Note that the uncertainty due to non-perturbative contribution is tiny compared to the scale error as we will show in the following discussions.The errors from ∆Q * and the UHO contribution decrease rapidly as the order increases.This has been shown more clearly in Figure 5.These two equations show that improved convergence can be obtained after using the PMCs scale-setting procedure.Table 5 shows that there are residual scale dependence for every order of the PMCs prediction, which, however, are markedly smaller than the conventional scale dependence of the conventional prediction.Note that with fixed √ s/Λ = 8, the scale variation in √ s/2 < µ r < 2 √ s leads to ∼ 10.1% variation for the QCD correction of R (4) uds (s), and ∼ 0.8% variation for the whole four-loop prediction R (4) uds (s).The PMC series, which has smaller (residual) scale dependence and a good convergent behavior, can be treated as the intrinsic perturbative nature of the series.
It should be pointed out that there are extra "power-correction" O(1/Q 4 ) to the total cross section in e + e − annihilation, which accounts for contributions that are fundamentally non-perturbative [64,65].These are introduced via non-vanishing vacuum expectation values originating from quark and gluon condensation.The non-perturbative addition to the Adler function has been calculated [64,66], where the non-perturbative operators are the gluon condensate, (α s /π)GG , and the quark condensates, m f qf q f .The latter obey approximately the partially conserved axial-vector current relations [67][68][69], where f π = (92.07± 0.85) MeV [3] is the pion decay constant.The complete dimension D = 6 and D = 8 operator are parameterized phenomenologically using the vacuum expectation values O 6 and O 8 , respectively.Note that in zeroth order α s , i.e. neglecting running quark masses, non-perturbative dimensions do not contribute to the integral, R(s) = Thus in the formula presented in Eq. (3.6)only the gluon and quark condensates contribute to R(s) via the ln s-dependence of the terms in first order α s .The gluon condensate cannot be fixed theoretically.There exist experimental determinations using finite-energy sum rule techniques: a fit using the τ vector plus axial-vector hadronic width and spectral moments yields, (α s /π)GG = (0.001 ± 0.015)GeV 4 [70], a moment analysis using cc resonances results in, (α s /π)GG = (0.017 ± 0.004)GeV 4 [71], an estimation on e + e − data gives the value of (α s /π)GG = (0.044 +0.004 −0.021 ) GeV 4 [72], a later fit on e + e − data yields, (α s /π)GG = (0.037 ± 0.019) GeV 4 [66].Due to the non-perturbative parameter (α s /π)GG fitted by different works are very different, the non-perturbative contribution will not be directly added to R uds (s), but just provided as a theoretical error, ±|R uds,NP | MAX , where the subscript "MAX" means the maximum of the absolute value |R uds,NP | when (α s /π)GG varying between the upper and lower bounds for all mentioned values, i.e., (α s /π)GG ∈ [−0.014, 0.056].Using these settings, we thus obtain a conservative estimate of the non-perturbative contribution for R uds (s), e.g., if taking the input parameter α s (M 2 Z ) = 0.1193 ± 0.0028 [73], yielding Λ = 352 +44 −41 MeV, we obtain ∆R uds ( √ s = 5Λ = 1.76GeV)|NP = ±0.0006,∆R uds ( √ s = 8Λ = 2.816GeV)| NP = ±0.00005.Remember that the non-perturbative contribution decreases rapidly as the centre-of-mass √ s increases, since it is proportional to 1/s 2 .There are thus total three theoretical errors for R uds (s) in our calculation: ∆R uds | ∆Q * , ∆R uds | UHO , and ∆R uds | NP .The total theoretical uncertainty is then determined by quadratically adding all the mentioned three errors.
It should be mentioned that our PMC calculation is based on the expression (2.3), which is obtained in the fixed-order perturbation theory (FOPT).The PMC procedure provides a resummation of all known higher order {β i } terms, thus a further improvement for the perturbation calculation of R(s).Another popular method to calculate R(s) is to evaluate numerically the contour-integral, R(s) , known as the contour-improved perturbation theory (CIPT) [74,75].The numerical solution of the contour-integral in CIPT involves the complete (known) RGE and provides thus a resummation of all known higher order logarithmic terms.The CIPT thus can also be applied to further improving the perturbation theory for R(s).The CIPT has also been widely used for the extraction of α s from τ lepton decays, e.g., Refs.[18,[76][77][78][79][80][81][82], and the extraction of α s from e + e − → hadron data, see, e.g., Ref. [83].The PMC calculation will also be used for the extraction of α s from R(s) data in next subsection.Further comparative investigation for the extraction of α s is left to future work.

Determination of α s
We adopt the PMC prediction (3.1) as the input to fit the R uds data in the energy range 1.84 GeV ∼ 3.72 GeV measured by KEDR Collaboration [11].All the data summarized in Table 16 of Ref. [11] are not independent but rather have point-by-point correlated effects, then the least squares (LS) estimators are determined by the minimum of χ 2 function [3], where e = (R exp.uds (s 1 ), R exp.uds (s 2 ), • • • , R exp.uds (s N )) is the column vector composed of N = 22 experimental data, and t = (R the.uds (s 1 ), R the.uds (s 2 ), • • • , R the.uds (s N )) is the corresponding column vector composed of theoretical predictions.The superscript T denotes the transpose.V −1 is the inverse covariance matrix which is derived from statistical errors and systematic uncertainties taking into account the correlation matrix presented in Table 18 of Ref. [11].The experimental uncertainty of the fitted parameter Λ is determined by requiring [3] χ 2 (Λ) = χ 2 min + 1. (3.9) The fitting results are presented in Table 6, where the first error is the experimental uncertainty.The second, third and fourth errors in 3rd column represent contributions from the residual scale dependence ∆Q * , the UHO of the pQCD approximant R ( ) uds | PMCs , and the non-perturbative power correction respectively.The second errors in 4th and 5th columns are determined by quadratically adding all the above mentioned three components, which represent the total theoretical uncertainty.All the components of the theoretical uncertainty are also determined by the minimum of χ 2 function (3.8), but calculating the theoretical predictions R the.uds (s i ) (= R ( ) uds (s i )) in the vector t = (R the.uds (s 1 ), R the.uds (s 2 ), • • • , R the.uds (s N )) by adding the corresponding theoretical errors, respectively.When calculating R the.uds (s i ) (= R ( ) uds (s i )) by taking various QCD corrections ( = 2, 3, 4) to fit the data, the running of the QCD coupling is changed accordingly, i.e, R MS and finally extract α s (M 2 Z ), as suggested by Ref. [84].We also present the value of χ 2 min , which represents the level of agreement between the measurements and the fitted function, and can be used for assessing the goodness-offit.For the 4-loop pQCD correction R uds (s)| PMCs with and the KEDR data is presented in Figure 6.The resultant α s (M 2 Z ) = 0.1227 +0.0117+0.0016−0.0132−0.0016 is consistent with the world average α s (M 2 Z ) = 0.1179 ± 0.0009 [3].Theoretical uncertainty is ∼ 1.3%, and is negligible compared to the experimental one (∼ 10%).Thus the accurate theoretical prediction (3.1) for R uds (s) allows to extract α s with high precision at the future Tau-Charm facility, such as the Super Tau-Charm Facility in China [85][86][87].It is necessary to emphasize that when using the conventional prediction to fit R uds data, the fitted Λ should satisfy the self-consistent requirement √ s/Λ > exp(π/2) 4.81.If using the conventional 4-loop prediction (3.3) to fit the KEDR R uds data [11], one may obtain an abnormal χ 2 (Λ) curve, see Figure 7, and an exaggerated Λ, thus all 16 data points below 3.3 GeV shall violate this constraint.Such violation can be highly improved if using the 4-loop PMCs prediction (3.1) to fit the data, where only the first 2 data points below 2 GeV violate this constraint.Note that the 2-loop and 3-loop fit will not violate this constraint, but will have larger theoretical errors.

Contributions to
Z ) Hadronic vacuum polarization (HVP) is not only a critical part of the Standard Model (SM) prediction for the anomalous magnetic moment of the muon, a µ = (g − 2) µ /2, but also a crucial ingredient for global fits to electroweak (EW) precision observables due to its contribution to the running of the fine-structure constant encoded in ∆α had (q 2 ).Traditionally, the leading order HVP contribution to a µ can be determined via the dispersion relation [88,89] where s th = m 2 π , α(0) is the fine-structure constant in the Thomson limit, the kernel function K(s) can be expressed analytically [88,89].
The running (scale-dependent) QED coupling, α(q 2 ) is determined via, α(q 2 ) −1 = α(0) −1 1 − ∆α had (q 2 ) − ∆α lep (q 2 ) , (3.11) where the contributions to the running are separated into hadronic (had) and leptonic (lep) components.The effective QED coupling at the Z boson mass, α(M 2 Z ), is the least precisely known of the three fundamental electro-weak (EW) parameters of the SM (the Fermi constant G F , M Z and α(M 2 Z )), and its uncertainty from hadronic contributions hinders the accuracy of EW precision fits.The hadronic contributions to α(q 2 ) are determined from the dispersion relation where P indicates the principal value of the integral.Using the PMC prediction (3.1), we evaluate the contribution of R uds (s) to a HVP,LO µ in energy range 1.8 ∼ 3.7 GeV, and present it as a function of Λ in Fig. 8, where the band represents the theoretical uncertainty, including contributions from the residual scale dependence and the UHO of the pQCD approximant for R uds (s) (3.1).Such theoretical uncertainty is ∼ 0.02% at Λ = 0.1 GeV and increases to ∼ 0.5% at Λ = 0.7 GeV.As for numerical results, if taking the same input α s (M 2 Z ) = 0.1182 ± 0.0012 as the KNT18 [90], we obtain a HVP,LO µ [1.841 ≤ √ s ≤ 2.00GeV] × 10 10 = 6.38 ± 0.02, where the total uncertainty includes effects from the α s uncertainty, the residual scale dependence, the UHO contribution, and the non-perturbative contribution.This result is in good agreement with the one reported by KNT18 [90], a HVP,LO µ [1.841 ≤ √ s ≤ 2.00GeV]×10 10 = 6.38±0.11,but with a decreased error, whose error is dominated by the variation of the renormalization scale µ r in the range √ s/2 < µ r < 2 √ s.When taking the same input as the DHMZ19 [91], i.e. α s (M 2 Z ) = 0.1193 ± 0.0028 from the fit to Z precision data [73], we obtain, where the error is obtained by quadratically adding the uncertainties from the α s uncertainty (±0.14), the residual scale dependence (±0.04), the UHO contribution (±0.01), and the non-perturbative contribution (±0.00).As for the running electromagnetic coupling at M Z , our prediction for the hadronic contribution from 1.8 ∼ 3.7 GeV range to the running of α(M whose uncertainty is dominated by the α s uncertainty (±0.10).The residual scale dependence (±0.02), the UHO contribution (±0.00), and the non-perturbative contribution (±0.00) are quite small.Our present predictions (3.13) and (3.14) are in good agreement with the results reported by DHMZ19 [91] but with decreased errors.

Summary
The hadronic e + e − annihilation rate R(s) is one of the most precise and theoretically safe observables involving strong interactions.The PMC provides a systematic method for solving the conventional renormalization scheme-and-scale ambiguities, and its PMC scale reflects the virtuality of the underlying QCD subprocess.By applying the PMCs, we have shown that a reliable and self-consistent analysis for R uds (s) can be achieved.Our new calculation for R uds (s) leads to a scale-invariant prediction, a significant stabilization of the perturbative series, and a reduction of theoretical uncertainty.It thus can provide a reliable and competitive determination for the QCD running coupling at future high-precision measurement on R uds (s), and will help to improve the accuracy of the SM predictions for the muon magnetic anomaly a µ as well as the QED coupling α(M 2 Z ).

Figure 4 .
Figure 4.The PMCs prediction R ( ) uds (s)| PMCs ( = 2, 3, 4) as a function of √ s/Λ.The red dotted, green dashed and blue solid curves are for the predictions at 2-loop, 3-loop, and 4-loop, respectively.For each curve, its error band represents the uncertainty from the residual scale dependence.
= 8.The errors of PMCs predictions are for the residual scale dependence, ∆Q * .The central values of conventional predictions (Conv.)are for µ r = √ s, and the errors are for ∆µ r .present the values of individual-order QCD correction terms for the fourloop predictions R (4) uds (s) with fixed √ s/Λ = 8 in
to -loop α s -running.For the computation of α s (M 2 Z ) based on Λ = Λ

Figure 6 .
Figure 6.Comparison of the fitted curve and the KEDR [11] data with statistical and systematic errors added in quadrature.The black thin line is for R (4) uds (s)| PMCs and its band is plotted by taking the fitted value Λ = 406 +207 −186 MeV.

Figure 8 .
Figure 8. a HVP,LO µ [1.8 ≤ √ s ≤ 3.7GeV] as a function of Λ.The band shows the total theoretical uncertainty, including effects from the residual scale dependence ∆Q * , the UHO contribution, and the non-perturbative contribution.

Table 1 .
The predicted 95.5% credible intervals (CI) for the UHO coefficients of ln Q 2 * /s via the Bayesian-based approach.The exact values (EV) are presented as a comparison.

Table 4 .
R = 5, 8, 12, respectively.The first errors are for ∆µ r of conventional predictions (Conv.), or ∆Q * of PMCs predictions (PMCs), the second errors are for UHO of the pQCD approximant R

Table 5 .
Total and individual-order QCD corrections for R

Table 5 .
The relative importance among

Table 6 .
[11]fitted Λ from R uds data below the D D threshold measured by KEDR collaboration[11].Results for R the.uds = R PMCs by taking various QCD corrections ( = 2, 3, 4) are presented.The first error is experimental, the second, third, and fourth errors of the 3rd column are from the residual scale dependence ∆Q * , UHO, and the non-perturbative contribution, respectively.The second errors of the 4th and 5th columns are total theoretical uncertainties, which include all the mentioned three theoretical errors added in quadrature.