Precise determination of αs from relativistic quarkonium sum rules

We determine the strong coupling αs(mZ ) from dimensionless ratios of roots of moments of the charm- and bottom-quark vector and charm pseudo-scalar correlators, dubbed RqX,n≡MqX,n1n/MqX,n+11n+1,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {R}_q^{X,n}\equiv \left({M}_q^{X,n}\right)\frac{1}{n}/{\left({M}_q^{X,n+1}\right)}^{\frac{1}{n+1}}, $$\end{document} with X = V, P , as well as from the 0-th moment of the charm pseudo-scalar correlator, McP,0.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {M}_c^{P,0}. $$\end{document} In the quantities we use, the mass dependence is very weak, entering only logarithmically, starting at 𝒪αs2.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\alpha}_s^2\right). $$\end{document} We carefully study all sources of uncertainties, paying special attention to truncation errors, and making sure that order-by-order convergence is maintained by our choice of renormalization scale variation. In the computation of the experimental uncertainty for the moment ratios, the correlations among individual moments are properly taken into account. Additionally, in the perturbative contributions to experimental vector-current moments, αs(mZ) is kept as a free parameter such that our extraction of the strong coupling is unbiased and based only on experimental data. The most precise extraction of αs from vector correlators comes from the ratio of the charm-quark moments RcV,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {R}_c^{V,2} $$\end{document} and reads αs(mZ) = 0.1168±0.0019, as we have recently discussed in a companion letter. From bottom moments, using the ratio RcV,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {R}_c^{V,2} $$\end{document}, we find αs (mZ) = 0.1186±0.0048. Our results from the lattice pseudo-scalar charm correlator agree with the central values of previous determinations, but have larger uncertainties due to our more conservative study of the perturbative error. Averaging the results obtained from various lattice inputs for the n = 0 moment we find αs(mZ) = 0.1177±0.0020. Combining experimental and lattice information on charm correlators into a single fit we obtain αs(mZ) = 0.1170±0.0014, which is the main result of this article.


Introduction
The strong coupling α s is the central quantity governing quantum chromodynamics (QCD). It is a key parameter to all observables computed in perturbation theory relevant for facilities such as the LHC or future e + e − colliders, which have an extensive program for determining top-quark and Higgs-boson properties such as their masses and couplings. It also plays a central role in flavor physics and in the determination of the masses of charmonium and bottomonium bound states. This parameter is also crucial for searches of physics beyond the Standard Model since it largely determines the size of the associated background. For a review on recent progress see e.g. refs. [1,2].

JHEP03(2020)094
A powerful method to determine parameters related to the strong interactions such as quark masses and α s are QCD sum rules based on weighted integrals of the total hadronic cross section R qq (with q = c, b) R qq (s) = 3s 4πα 2 σ e + e − → qq +X (s) σ e + e − → qq +X (s) σ e + e − → µ + µ − (s) . (1.1) Particularly important for our work are the inverse moments, M V,n q , of R qq (s) defined as M V,n q = ds s n+1 R qq (s) . (1.2) Using analyticity and unitarity, these can be related to the coefficients of the Taylor expansion of the quark vector-current correlator around s = 0, which can be computed rigorously in perturbative QCD (pQCD) for n not too large. A shortcoming of using moments M V,n q is that, while the integration in eq. (1.2) over the normalized cross section extends all the way to infinity, experimental data are limited to a finite energy range. If the energy of the last measured cross section is sufficiently large, one can safely use the theoretical prediction for the R-ratio in perturbation theory as a substitute (the region is sometimes referred to as the continuum), applying some penalty to reduce the model dependence. For the charm cross section the data above threshold spans a wide range of energies such that even for n = 1 the computed moment is fairly insensitive to how the continuum is treated [3]. On the other hand, bottom moments with low values of n do depend strongly on the continuum such that M V,1 b cannot be used for any competitive determination of the bottom-quark mass [4,5] -a situation that could change if data at larger energies became available. Here, since we are interested in a precise extraction of α s , the continuum contribution must be treated carefully, in a way that avoids any possible contamination of the extracted values.
Theoretically, the moments M X,n q are governed by the typical scale m q /n Λ QCD . This is easy to understand since large values of n have more weight in the narrow resonances such that a non-relativistic treatment becomes necessary. For small values of n one can compute the moments in perturbative QCD supplemented by non-perturbative power corrections parametrized in terms of local condensates. This framework is known as the operator product expansion (OPE) [6,7]. It turns out that the perturbative term overly dominates the series (even more so for the bottom quark) and the leading (gluon) condensate is introduced mainly as an estimate of the size of non-perturbative corrections. This method goes under the name of relativistic quarkonium sum rules.
An interesting alternative which does not suffer from problems related to the continuum are moments of the pseudo-scalar quark-current correlator, which can be accurately computed in lattice QCD [8] -although, so far, precise simulations exist only for the charm quark. Interestingly, the 0-th moment of this correlator is physical, 1 and quite insensitive to the charm-quark mass, which makes it an ideal candidate to determine α s .

JHEP03(2020)094
On the other hand, it has been shown that the perturbative series of the pseudo-scalar moments (at least for n > 0) displays a quite poor convergence [5].
A lot of progress has been made in the lattice community for determining QCD parameters from the pseudo-scalar correlator since the pioneering work of ref. [8], in which the charm-quark mass and the strong coupling were extracted (the former with high accuracy). Focusing on α s , the follow-up paper by HPQCD [9] already claimed half-percent accuracy at the Z-boson mass with a value very close to the world average, while refs. [10] and [11] have somewhat smaller central values and slightly larger uncertainties (0.7 % and 1%, respectively). The latter references introduced an interesting alternative strategy to determine α s : the use of ratios of moments for which the mass dependence largely cancels; a strategy that we extend and exploit here also for the vector-current moments.
So far, nearly all lattice analyses with sub-percent accuracy on α s (m Z ) estimate the perturbative uncertainties in essentially the same way: fixing the renormalization scale to some default value and adding an estimate for the first unknown perturbative coefficient. The truncation error is estimated by either varying the size of the guessed term or comparing the values of the strong-coupling constant obtained including or not the next term (an exception to this paradigm is the analysis of ref. [12]). However, as argued in refs. [3,5], the renormalization scale of m q plays an important role in realistic estimates of perturbative uncertainties. In ref. [12], which only uses ratios of moments, perturbative uncertainties are estimated performing renormalization scale variation in a manner analogous to refs. [3,5]. In particular, renormalization scales for α s and the heavy-quark mass are floated independently within a certain range. Their result is compatible with the world average, but the quoted uncertainties are more conservative, achieving only 2.2% accuracy. This shows that, while it can probably be argued that the computation of the moments on the lattice is fairly under control, one could question some of the estimates of theory incertitudes found in the literature.
In this article we extract the strong coupling with n f = 4 and n f = 5 from QCD sum rules analyzing dimensionless mass-insensitive ratios of the moments M V,n q and M P,n q , denoted R X,n q and defined in eq. (2.6). This type of strategy was originally introduced for pseudo-scalar moments computed on the lattice [10,11] and, to the best of our knowledge, was applied to vector-current correlators, whose moments can be obtained using real experimental data on R qq (s), for the first time in a companion paper [13]. Here, we provide further details on the analysis of ref. [13], extend the method to bottom vector-current moments, and, finally, apply the same strategy to the lattice-determined pseudo-scalar moment ratios of the charm current.
On the experimental side, we compute the ratios R V,n q using data on narrow resonances and the total hadronic cross section above the charm and bottom open thresholds, supplemented with perturbation theory in the regions where no measurements exist (and to subtract the uds background from charmonium data). Since our goal is to perform a rigorous extraction of α s , the value of the coupling used as input in the perturbative expressions must be kept as a free parameter, such that α s is extracted from the experimental data without any bias. On the theoretical side, we carefully examine perturbative uncertainties and conclude that the renormalization scales of α s and the MS quark mass have to be var-

JHEP03(2020)094
ied independently within a certain range to avoid theoretical biases. Using this approach, we also perform a reanalysis of lattice data on the ratios of moments R P,n c and the M P,0 c regular moment for the pseudo-scalar charm correlator. The method provides a reliable extraction of α s and, in particular for the vector-current charm-quark ratios, R V,n c , the values obtained for α (n f =5) s (m Z ) are competitive given the accuracy of current experimental data. Finally, results from lattice QCD and the charm vector moment ratios can be unambiguously combined to obtain an even more precise strong coupling determination. This paper is organized as follows: in section 2 we provide an overview of the theoretical ingredients that enter our analysis, both perturbative (section 2.1) and non-perturbative (section 2.2); we review the experimental and lattice input which is used to determine α s in section 3; a detailed study of perturbative uncertainties is given in section 4, while our results are presented in section 5 and compared to previous determinations in section 5.3; finally, our conclusions are contained in section 6.

Theoretical input
In this section we discuss the theoretical description of inverse moments of the vector and pseudo-scalar quark-currents, as well as the ratios formed from these that we exploit in the present work. The moments of eq. (1.2) can be related, using analyticity and unitarity, to the Taylor coefficients of the expansion of Π V q at s = 0 as [6,7] with √ s = p 2 , the e + e − invariant mass, Q q the quark electric charge, q = c, b, and is the quark vector current.
Using the notation of ref. [5], we define the pseudo-scalar quark-current correlator as with j P q (x) = 2 m q iq(x)γ 5 q(x); here we will only consider pseudo-scalar moments of the charm-quark current (q = c). The additional mass factor in the pseudo-scalar current (as compared to the vector case) makes it formally scheme and scale independent. Moments analogous to those of eq. (2.1) can be defined as where we use the combination first introduced in ref. [5]

JHEP03(2020)094
The theoretical quantities that will be used in this article to determine α s are mass insensitive and dimensionless. For the case of the pseudo-scalar correlator one can use the 0-th moment, which has mass dimension zero by itself, and depends on the quark mass only logarithmically starting at O(α 2 s ). This moment is an observable, in the sense that it does not need an ultraviolet subtraction to become finite, being formally renormalizationscale and scheme independent (although it still retains a residual µ dependence at any finite order in perturbation theory). The 0-th moment of the vector correlator, on the contrary, cannot be related to any experimentally measurable quantity. It is related to the subtraction that renders the sum rule finite, and is therefore scheme and scale dependent.
One can, however, work with ratios of roots of the n > 0 moments, such that the mass dependence almost completely disappears. The quantities we are interested in are the ratios of consecutive roots of moments. Specifically, we define the following massinsensitive quantities where X = V, P refers to vector and pseudo-scalar correlators, respectively. This type of ratio of moments was originally introduced for the pseudo-scalar correlator [10,11]; here we extend their use to the vector-current as well. They are the central objects of our analyses.
We write the perturbative vacuum polarization function for vector and pseudo-scalar currents expanded around s = 0 as (2.7) To have a common notation for both currents we use Π P q (q 2 ) = P q (q 2 ), where P q is the twice-subtracted pseudo-scalar correlator defined in eq. (2.5). 3 In full generality, different renormalization scales can be employed for the mass and the coupling in the perturbative expansion of the moments. We denote those scales µ m and µ α , respectively. As shown in refs. [3,5], in order to properly assess the size of perturbative uncertainties, it is important to vary them independently. However, for the 2 In ref. [20] the three-loop vector correlator has been obtained numerically for any value of s/m 2 to arbitrary precision. 3 To simplify our nation, here and in what follows, we do not write explicitly the dependence on the number of flavors n f since it can be deduced from the context.  time being, it is sufficient to set both scales to the quark mass µ m = µ α = m q , employing the shorthand notation m q ≡ m q (m q ). With this choice, the logarithms are resummed and the perturbative expansion of the moments in powers of α s takes the following simple form The ratios of moments we are interested in, defined in eq. (2.6), are dimensionless and mass insensitive. Their perturbative expansion in powers of α s with the choice µ m = µ α = m q can be written aŝ It is convenient to organize the computation of the coefficients r X,n i taking first the logarithm of the moment, expanded in powers of α s as log (2 m) 2nM X,n q = log c X,n with coefficients that obey the following recursive relation in terms of the original ones c X,n 0 a X,n i+1 = c X,n i+1 − The logarithm of the ratio is now trivially expressed as The last step to obtain the fixed-order series for the ratios of moments is expanding eq. (2.12). This can be done using the following computer-friendly recursion relation (2.13)   Table 3. Perturbative coefficients M P,0 i and r P,n i for the 0-th moment (second row) and the ratios of moments (rest of rows) for the charm pseudo-scalar current correlator. We show only terms which do not involve logarithms.
Since we are interested in assessing the size of perturbative uncertainties through renormalization scale variation, we need to express eq. (2.6) in terms of α s (µ α ) and m q (µ m ). In a first step to compute the associated logarithms, we set both scales equal and definê , (2.14) with r X,n i ≡ r X,n i,0 . Imposing that eq. (2.14) does not depend on µ m one can obtain r X,n i,j in terms of r X,n i defined in eq. (2.9), and the QCD beta function and MS mass anomalousdimension coefficients, which are defined, respectively, by They are known to five loop accuracy [30][31][32][33][34][35][36][37], and collected in table 4. In terms of those we find An important feature of this perturbative expansion is that the first non-zero coefficient of a logarithmic term is r X,n 2,1 and, therefore, the dependence on m q starts only at O(α 2 s ). This weak logarithmic dependence on m q makes our results rather insensitive to the quark mass.

JHEP03(2020)094
The most general formula has in general µ α = µ m , depends on two kind of logarithms, and contains two nested sums: with r X,n i,j,0 ≡ r X,n i,j and where we made explicit that the first logarithm appears only at order α 2 s . We can easily obtain the r X,n i,j,k in terms of r X,n i,j and the QCD beta function coefficients imposing that eq. (2.17) is independent of µ α . This results in the recursive formula first introduced in [38] For the 0-th moment of the pseudo-scalar correlator M P,0 c , eqs. (2.14) to (2.18) still apply with the obvious replacement r X → c P . In some occasions it will be convenient to abuse notation defining R P,0 c ≡ M P,0 c . The numerical values for the coefficients of ratios of vector moments are collected in tables 1 (charm) and 2 (bottom). For the charm pseudo-scalar correlator, the perturbative coefficients are given in table 3.
The total α s correction at order O(α 3 s ) to the first three charm-quark vector-current ratios is of about 12.5% for R V,1 c , 7.2% for R V,2 c , and 5.2% for R V,3 c . The bottom-quark vector-current ratios are less sensitive to α s corrections, which turn out to be almost a factor of two smaller than in the charm-quark counterparts. The first ratio, R V,1 b , receives a correction of about 7.7%, while for the next two ratios the correction is 4.1% and 2.8%, respectively. Precise extractions of α s from bottom quark ratios require, therefore, smaller experimental errors in order to overcome the smallness of pQCD corrections. The 0-th charm pseudo-scalar moment is again quite sensitive to pQCD, with a total correction of about 29%. The first two ratios of pseudo-scalar moments receive large α s corrections as well: 24% for R P,2 c and 11% for R P,2 c . This higher sensitivity is welcome in the sense that less precision is required from the lattice computations. The price to pay, however, is that the pseudo-scalar moments display a bad perturbative convergence which leads to larger errors from the truncation of the perturbative series.
All the formulas and recursive relations in this section have been implemented into a numerical python [39] code. We have also written an independent Mathematica [40] program which is based on a direct derivation of the formulas using built-in Mathematica functions. Our codes agree within 15 decimal places. This completes the description of the perturbative contribution to the moments and ratios of moments.

Non-perturbative corrections
The perturbative contribution presented in the previous section will be corrected for nonperturbative effects including the most important sub-leading contribution from the OPE, which is the gluon condensate. Since this matrix element has mass-dimension 4, its contribution is suppressed by four powers of the heavy-quark mass [41,42]. The renormalizationgroup invariant (RGI) scheme for the gluon condensate [43] will be used in our analyses.  Table 4. QCD beta function and MS mass anomalous dimension coefficients for n f = 4 (second and fifth rows) and n f = 5 (third and last rows), up to five loops.
For the RGI gluon condensate we take the value [44] α s π G 2 As argued in ref. [45], it is convenient to express the gluon condensate Wilson coefficient in terms of the pole mass to stabilize the correction for large values of n. To obtain a numerical value for the pole mass m (pole) q we follow ref. [3] and use the one-loop conversion: Therefore, in practice the pole mass depends both on µ m and µ α . For the purpose of this section (that is, to obtain the non-perturbative corrections to R X,n q ), we consider the series for the n-th moment only at O(α s ) for both perturbative and non-perturbative terms, which we write schematically as where, for notation simplicity, we do not explicitly show the dependence of C X,n 1 on µ m and m q (µ m ). The gluon-condensate Wilson coefficients, g X,n i , for current-current correlators are known to O(α s ) [46]. The numerical values for the charm vector correlator can be found in table 5 of ref. [3], while table 6 of ref. [5] shows the values for the bottom vector and the charm pseudo-scalar moments. Next we take the n-th root and expand in α s and the gluon condensate up to linear terms, which gives n . (4) [α/α(M )] 2 0.957785 0.95554 Table 5. Masses and electronic widths [47] of the narrow charmonium resonances and effective electromagnetic coupling [48]. α = 1/137.035999084(51) is the fine structure constant, while α(M ) stands for the pole-subtracted effective electromagnetic coupling at the scale M .
From this we can take the ratio of two consecutive moments. We define a n = (c n,X 0 ) 1 n , b n = a n n C X,n 1 , where, again, we refrain from explicitly displaying the dependence of a n , b n , c n , and d n on the current, µ m , m q (µ m ) and m (pole) q . In terms of those, the gluon condensate correction to the ratios can be finally written as .
For the n = 0 pseudo-scalar moment M P,0 c one directly uses the gluon condensate correction shown in eq. (2.21) with g P,0 0 = 4/15 and g P,1 0 = 1.4086. The relations derived in this section have been implemented into our python code, while a direct computation is included into our independent Mathematica program. This concludes the presentation of the theoretical input. We give further details of our numerical codes in section 4.

Experimental and lattice input
In this section we briefly discuss how to obtain the experimental values for the moments that go into our analyses. For the pseudo-scalar correlator there is, of course, no experimental information available, but several lattice collaborations have obtained numerical values for low-n moments, as well as ratios of moments, through numerical Monte Carlo (MC) simulations.

Charm vector correlator
Regular moments of charm-tagged cross section where already worked out in ref. [3], which used data on narrow charmonium resonances as given in [49], less precise than what is available today. In this section we only streamline the procedure, and update our results -10 -

JHEP03(2020)094
using the values for the charmonium electronic widths provided in ref. [47], which are collected in table 5. 4 There are three main contributions to the moments: narrow resonances (which appear below threshold); threshold (with data collected by different experimental collaborations [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66] up to energies of 10.538 GeV); and continuum, from 10.538 GeV to infinity, with no experimental data available. In ref. [3] a method to recombine all available experimental data into a single dataset was employed. It was based on the algorithm used in ref. [67] to reconstruct the total hadronic cross section below the charm threshold, in the context of the computation of the vacuum polarization function contribution to the g − 2 of the muon. 5 The method was generalized to subtract the non-charm background, which uses α s -dependent theoretical predictions for the light-quark and secondary charm production cross sections. In the original implementation of [3] this pQCD prediction was re-weighted by a constant factor determined from a fit to data below the open charm threshold. At that time, the most precise data in that region was provided by BES [52,53,[55][56][57] inclusive measurements, which are much higher than theoretical pQCD predictions and data from exclusive measurements. Newer experimental measurements by KEDR [69][70][71] are significantly lower and compatible with pQCD. Hence, in the present analysis of charm data, we directly use pQCD without any additional normalization. Using this subtraction ansatz, the perturbative QCD prediction for the charm-tagged cross section R cc at energies E ∼ 10 GeV is in full agreement with experimental measurements.
We have re-written our old Mathematica program into a python code, which uses iminuit [72], the Minuit [73] python implementation. Our updated code exactly reproduces the results in ref. [3] if the same input is used, but it yields more precise results for the moments when the experimental values for the narrow-resonance electronic widths are updated. More importantly, it also computes the correlation matrix among the moments, which is essential to correctly determine the uncertainties for the ratios of moments.
Finally, since we aim at a precise determination of α s and the moments depend on this parameter through the non-cc background 6 and the continuum, we should obtain the ratios of moments R V,n c as a function of α s . This allows for an unbiased extraction of the coupling with the continuum and background contributions determined self-consistently. We have then computed the moments for multiple values of α s , finding a remarkably linear dependence on the coupling which allows for an accurate (and simple) parametrization of R V,n c in terms of α (n f =5) s (m Z ). We find that the dependence of the moments with α s is monotonically decreasing, due to the higher weight of the non-charm background 4 These correspond to the values used in ref. [50] which is an update of a previous analysis on the determination of the charm-quark mass. We do not adopt the current PDG average [51] (2019 update), that quotes slightly different results for Γee: 5.53(1) keV and 2.33(4) keV, which are fully compatible, with slightly larger uncertainties for the J/ψ. Using these values our results for R V,n=1,2,3 c would decrease by 0.1%, 0.06% and 0.04%, shifts which are 10, 4 and 3 times smaller than our uncertainties on those. The effect on the fitted value of αs(mZ ) would be, in the worse case, a 0.2% shift downwards, 10 times smaller than the total uncertainty. Furthermore, we have not updated the value of the charmonium masses since the uncertainty is overly dominated by the electronic widths. 5 An alternative treatment of the contribution above the charm threshold can be found in ref. [68]. 6 The small O(α 3 s ) singlet contribution [74], as argued in ref. [48], is very small and can be safely neglected for both charm and bottom analyses.  subtraction as compared to the continuum contribution. The moments' uncertainty is dominated by data and found to be α s independent. We quote the results obtained for the ratios R V,n c as a function of α 7 These results were reported for the first time in ref. [13]. Since the uncertainties among the various moments are highly positively correlated, the ratios turn out to be more precise than the individual moments. While the relative precision for the first 4 moments is roughly constant and around 1%, the uncertainties for the first 3 ratios rapidly decrease as n grows giving 0.98%, 0.22% and 0.104%, respectively. This is partially caused by the fact that the narrow-resonance contribution (with very small errors) has a stronger weight for larger n. The value for the higher ratios seems to freeze and we find R V,n→∞ c → 1.

Bottom vector correlator
Regular moments of bottom-tagged cross section where discussed in detail e.g. in ref. [5]. In this case one has to combine the contribution from the first four narrow resonances with threshold data from BABAR [75], which has to be corrected for initial-state radiation and vacuum polarization effects. This unfolding of the data induces a correlation among the different data points, which in turns translates into a stronger correlation for the moments. We have translated our old Mathematica program that performs the QED corrections into a fast python code, which allows to take many more iterations in the unfolding procedure using very little CPU time, accurately reproducing the results of ref. [5]. BABAR data stops at 11.52 GeV, and some modeling becomes necessary at larger energies. The approach of ref. [5] was to interpolate the last experimental points with the pQCD prediction in a smooth way, assigning an energy-dependent systematic uncertainty to the model, linearly decreasing with the invariant squared mass s from 4% at Q = 11.52 GeV to 0.3% at Q = m Z . Here, we tune the dependence of the uncertainty with s in accordance with expectations from hadronization power corrections based on the operator product expansion, parametrized by the gluon condensate, which predicts a dependence of the type 1/s 2 . This results in a moderate reduction of the moments' uncertainty.
Our updated code also provides the correlation matrix among the moments, which is used to calculate the ratios' uncertainties. Finally, since there is some (small) α s dependence left in the moments through the continuum [ this includes the perturbative QCD prediction and an interpolation between pQCD and a linear fit to the (QED corrected) BABAR data for energies larger than 11.05 GeV ], we again evaluate the moments for 7 The updated results for individual moments M V,n c will be given elsewhere.   The partial cancellation of correlated uncertainties in the ratios of bottom moments is much larger than for charm. While regular moments with n < 5 are rather imprecise, with relative accuracy of 1.45%, 1.38%, 1.26%, and 1.20%, respectively, the first three ratios are below the percent accuracy: 0.55%, 0.23%, 0.12%. Finally, we also observe the behavior R V,n→∞ b → 1 in bottom moments, but in this case the limit is approached from below.

Charm pseudo-scalar correlator
Although the pseudo-scalar current is not accessible in experiments in the same way as the vector (i.e. there is no such thing as a "pseudo-scalar photon"), results for the associated moments can be obtained from simulations on the lattice. The experimental input is effectively passed to the simulation by tuning lattice parameters to a number of physical observables. The tuned lattice action (which is no longer modified or adjusted) is then used to perform the predictions for the moments. Lattice simulations have to overcome some difficulties, such as the continuum, infinite-volume, and physical mass extrapolations, which can translate into sizable uncertainties. The simulations are based on MC methods, and are therefore also limited by statistics. Other aspects to take into account when using lattice data is which type of action is used to compute the fermion determinant. According to ref. [8], moments of the pseudo-scalar current are not as afflicted by systematic uncertainties as the vector-current ones, and therefore might be used for precision analyses.
Lattice results for the pseudo-scalar correlator are provided in terms of the so-called reduced moments R k . They are constructed to have mass-dimension 0, and factor out the tree-level results such that their perturbative expression starts with a coefficient equal to 1. For n = 0 and n > 0 they are related to the notation in eq. (2.7) as follows The mass dimension of "regular" moments with n > 0 is obtained through powers of the η c mass. However, when taking ratios this dependence, as expected, completely drops, and -13 -

JHEP03(2020)094
we obtain the following relation: Numerical values for the 0-th moment of the pseudo-scalar correlator, as well as for the first two ratios, are given in table 7. We have transformed the results quoted by various lattice studies into our notation. We again observe that systematic lattice uncertainties in regular moments largely cancel when taking the ratios, and, even though only the first two have been computed, one can see that central values decrease when going from n = 1 to n = 2, becoming closer to 1.

Analysis of perturbative uncertainties
In this section we investigate the convergence properties of those perturbative series that will be used to determine the strong coupling. For that we will need to evolve α s and the MS heavy-quark masses with the renormalization scales µ α and µ m , respectively, using the corresponding renormalization group evolution (RGE) equations. Details on how this is implemented are provided in appendix A.
In order to thoroughly study the perturbative uncertainties, and following refs. [3,5], we use two independent renormalization scales, which we call µ α and µ m . The perturbative series we shall be dealing with are written in terms of α s (µ α ) only. It is important to have a single expansion parameter [ that is, one has to avoid having α s (µ m ) explicitly in the series ], such that the pole-mass related renormalon is properly canceled. Therefore, the dependence on µ m starts only at α 2 s , as powers of log(µ α /µ m ) and log[ µ m /m q (µ m )]. Hence, it is expected that the dependence on µ m is weaker than on µ α , which in turn might mean that double scale variation is not as crucial as in quark mass determinations. In any case, to be conservative, we adopt the same scale variation as in [3,5]: m q ≤ µ α , µ m ≤ µ max , with µ max = 4 (15) GeV for charm (bottom). For practical purposes we create a grid in the two renormalization scales, with 4000 and 3025 evenly distributed points for bottom and charm respectively, which correspond to a bin size of ∼ 0.05 GeV. We will explore how uncertainties change if other conventions, some of which less conservative, are adopted.
When performing scale variations, it is customary to avoid extreme cases where the series converges badly or contains large logarithms. Firstly, we performed an analysis of the convergence properties of the perturbative series for α s in the spirit of what was done in ref. [5], studying the convergence of each series for different values of µ α and µ m within our grids. In ref. [5] it was suggested that series with bad convergence properties could be discarded. However, in the present case we find a rather flat distribution for the parameter that measures the convergence of the series, in contrast to what was found in ref. [5] for quark-mass determinations. The detailed results of this analysis are given in appendix B. Instead, here we shall use a more standard criterion based on avoiding large logarithms, and will simply require that 1/ξ ≤ µ α /µ m ≤ ξ, being our canonical choice ξ = 2. The excluded regions are shown as faint gray areas in figure 1. We do not impose a similar veto on µ m /m q (µ m ) since the original variation range on µ m already implements -14 -JHEP03(2020)094 the usual small-log paradigm (here one cannot use values of µ m smaller than m q since then α s becomes large and endangers the convergence properties of the series). Furthermore, the analyses in appendix B also show that the bottom vector ratios are the most convergent, closely followed by the charm vector correlator. However, the series for the pseudo-scalar correlator are significantly less convergent than the other two cases studied, a behavior that was already found in ref. [5]. Therefore, determining α s from the charm and bottom vector correlators seems warranted, at least from the perspective of perturbative uncertainties.
We continue our exploration of perturbative uncertainties drawing contour plots that show the dependence of the α (n f =5) s (m Z ) extracted value on the renormalization scales. For this exercise we do not include the gluon condensate correction and use the experimental values quoted in table 6 for the vector correlator (using the world average value for the strong coupling), and the results of [10] shown in table 7 for the pseudo-scalar correlator. We also use the values m c = 1.28 GeV, and m b = 4.18 GeV for the quark masses. We analyze the values of α s as obtained from the series which include up to O(α 3 s ) terms. The results for the various currents and number of flavors are collected in the three rows of figure 1, where different columns correspond to different ratios (except in the last row, where the leftmost panel shows the result for the n = 0 pseudo-scalar moment) and gray shaded areas show excluded values for the choice ξ = 2. From the plots one can conclude that, in most cases, varying the scales in a correlated way in some limited ranges may lead to serious underestimates of perturbative uncertainties. In some cases, however, a variation keeping both scales equal seems to capture most of the spread in α s values, while in others keeping µ m = m q seems also sufficient. It seems that there is no unique 1D slice valid for all situations, and hence we conclude that an independent variation of scales is the safest procedure.
We turn now to a systematic study of the ξ dependence. Since at O(α s ) there is no µ m dependence, to estimate the perturbative uncertainty at this order we simply vary µ α between m q and µ max . At higher orders we once again vary m q ≤ µ α , µ m ≤ µ max , with the constraint 1/ξ ≤ µ α /µ m ≤ ξ. The perturbative uncertainty is then computed as half of the difference between the maximum and minimum values of α s in the constrained grid, while the central value is simply the average of those two. The standard choice to estimate perturbative uncertainties is to vary the argument of logarithms dividing/multiplying it by factor of 2, which corresponds to our canonical choice ξ = 2. Taking ξ = 1 corresponds to the correlated variation µ α = µ m , while very large values of ξ do not impose any constraint. We show the dependence of the central value and perturbative uncertainty on the value of ξ in figure 2. To carry out this study we do not fix the value of α s on the charm and bottom vector-correlator experimental moments, since setting it to the world average makes the uncertainty smaller and we want to assess the effect of ξ on our final uncertainty. For the bottom vector correlator, taking both scales equal yields uncertainties smaller than our canonical choice by factors of 1.4, 2.4, and 2.3, for the first three ratios, respectively. For the charm-vector correlator the difference grows dramatically, with uncertainties that are smaller by factors of 3.7, 2.5 and 2.1 for the first three ratios. In the pseudo-scalar case the uncertainty as a function of ξ is nearly the same for all quantities considered, and the correlated estimate is 2.4 times smaller than the default choice. The unconstrained error  estimate is at most 53% larger than that with ξ = 2 for all cases. Except for values of ξ close to 1, the central value grows as ξ increases, but the variation is below the percent in all cases. We finish this section by exploring the order-by-order convergence of the extracted values of α s . Again we ignore non-perturbative effects and fix the quark masses. We also assume experimental moments have no uncertainties, such that error bars shown in this section are purely of perturbative origin. Taking the default constraint ξ = 2 we obtain the results shown in figure 3. We see excellent convergence between the O(α 2 s ) and O(α 3 s ) determinations in all cases, while for ratios with n > 1 there is a slight tension between the O(α 1 s ) and O(α 3 s ) results. This is not cause for concern since the LO extraction does not yet depend on µ m and therefore should be regarded as a special case. A similar situation was found with the O(α 0 s ) determination of quark masses in e.g. ref. [5], which was independent of µ α . In that sense the O(α n s ) quark-mass determination should be thought of as being of the same order as the O(α n+1 s ) strong-coupling extraction.

Results
In this section we present the main results of our analysis: the determination of α  The associated uncertainties are very small and barely contribute to the final error, since the ratios of moments we use are rather insensitive to the quark mass. Non-perturvative corrections also contribute to the error budget, but their contribution is absolutely negligible in the case of the bottom-quark based determinations, and always subleading for the charm-quark ones. For charm-quark based determinations one could consider an additional source of uncertainty coming from matching the theories with n f = 4 and 5 active flavors at the scale µ b , which by default is taken to be m b . The choice of µ b inflicts a tiny uncertainty, which we estimate by considering µ b = 2 m b and m b /2. Running α s at 5 loops (and matching accordingly at 4 loops) it turns out to be negligibly small: 5 × 10 −6 . The uncertainty on the bottom mass also affects the matching relation, but the associated error is also insignificant: 1 × 10 −5 . These are much smaller than any other source and will no longer be mentioned.

α s from ratios of vector correlators
In this section we present results based on "real" experimental data, that is, α s extractions from ratios of vector-correlator moments, for n f = 4 and 5. For these analyses we use the  Table 8. α (n f =5) s (m Z ) determination from ratios of bottom and charm vector-correlator moments R V,n q , eq. (2.6). The first column specifies the flavor content of the current, the second column shows which ratio has been used, while the third gives the central value. Fourth to seventh provide the various components of the uncertainty: scale variation (σ pert ), experimental (σ exp ), quark mass (σ mq ), and gluon condensate (σ np ), which are added in quadrature in the last column (σ tot ). α s dependence of the experimental moments, given in table 6, solving the relevant equations consistently. The determinations from the charm correlator are shown graphically in figure 4(a). For comparison, the world average is shown as a faint gray band. All charmonium (and bottomonium) determinations are compatible among themselves and with the world average. Both extractions are quite robust, with rather stable central values, although the extraction from bottomonium yields somewhat larger central values than the extractions from charmonium sum rules. A detailed splitting of all sources of incertitude is given in table 8. For both quarkonium systems we observe that perturbative uncertainties grow with n (this behavior was already seen in figure 3), particularly for charmonium, with overall larger errors. Experimental uncertainties behave in the opposite way, decreasing as n grows. This is expected since larger values of n are dominated by the very precise narrow-resonance contribution. For charmonium, the larger experimental uncertainties discards the first ratio for precision extractions of α s . If the experimental error could be drastically reduced, n = 1 could however turn into a competitive measurement, since, from the theoretical point of view, it is quite clean. For both systems there is a compensation of the two effects such that the uncertainties for n = 2 and n = 3 are comparable. Since the ratios R V,2 q involves the moments M V,2 q and M V,3 q their perturbative expansion is expected to be better behaved than the ones for R V,3 q which brings the contribution from M V,4 q . Larger values of n are disfavored since non-relativistic duality violations could start playing a non-negligible role. This is in line with the smaller perturbative uncertainty for R V,2 q . Since the charm-quark based extractions have smaller errors and in the spirit of avoiding moments with large values of n we take the charm n = 2 result as the main outcome from the analysis of vector correlators:   [11], and Nakayama et al. [12]. n = 0 corresponds to R P,0 c ≡ M P,0 c , while n = 1, 2 stand for R P,n c . Perturbative uncertainties are estimated trimming the grids with the default parameter ξ = 2. The light gray area shows the current (2019) world average with its uncertainty [47].
being fully compatible with it: the central values differ by 0.6 σ. This value of α s was reported, for the first time, in ref. [13].
It is interesting to compare our final uncertainties, based on the conservative procedure of varying both scales independently, with more optimistic approaches often used in the literature. For example, if we had used a correlated scale variation with µ α = µ m the central value would decrease by 0.0001, but the perturbative uncertainty would decrease by a factor of 2.12 to (0.0007) pert , making up for a total error of only (0.0013) total . On the other hand, using a completely unconstrained grid, the central value and perturbative uncertainty grow to 0.1175 ± (0.0022) pert , a 50% increase in error, with a total uncertainty of (0.0025) total .

α s from lattice data
We turn now our attention to the pseudo-scalar Green's function, for which "experimental" data is obtained from lattice MC simulations. A number of collaborations provide results for the same quantities, and we analyze all of those with the same theoretical expressions and identical treatment of perturbative uncertainties. The results are shown graphically in figure 4(b), using different colors for the various lattice determinations. From left to right these are: Allison et al. [8] and McNeile et al. [9] (HPQCD collaboration, MILC configurations, HISQ action for light quarks and c-quark propagator, only M P,0 c ); Maezawa et al. [10] and Petreczky et al. [11] (HotQCD configurations, 2 + 1 flavors plus valence cquarks treated with the HISQ action, M P,0 c and R P,n c with n = 1, 2); and Nakayama et al. [12] (JLQCD collaboration, 2 + 1 flavors plus valence c-quarks treated with the Möbius -20 -
We observe that the total uncertainty is overly dominated by the truncation error, such that efforts to compute these quantities on the lattice more precisely are not warranted with our present knowledge of the perturbative series. Interestingly, perturbative errors for the pseudo-scalar correlator seem roughly independent of the moment being used. Even though all determinations are well compatible with the current world average, we observe that results from R P,2 c are slightly lower than those with M P,0 c and R P,1 c . This behavior is less pronounced for the JLQCD results, which could suggest that the shift is caused by lattice results (another argument in favor of this reasoning is that for the vector correlator results with larger n are higher). We also observe that M P,0 c -based extractions are higher if HPQCD results are employed, although still compatible. For R P,1 c -based determinations, the newest HotQCD result of ref. [11] is in very nice agreement with the JLQCD result, although the old HotQCD determination [10] is significantly larger. Interestingly, for determinations with R P,1 c the situation is the opposite, and there is excellent agreement between JLQCD and the old HotQCD results, being again the extraction of [11] somewhat lower. All in all, the various results are compatible and there seems to be a higher density of central values around the world average. As a representative value of α s from the pseudoscalar lattice correlator we average the various central values from the n = 0 moment, finding α We finish this section by comparing the perturbative uncertainties for vector and pseudo-scalar correlators with charm quarks. As expected from the analysis carried out in section 2.1, the uncertainties for the pseudo-scalar correlator are larger for the n = 1, 2 ratios of moments, but for n = 3 they become of the same order. This is in line with the findings of ref. [5] for regular moments, and again points to the fact that the total uncertainty will not go down with more precise lattice simulations. Without 5-loop results, possible ways of improving the accuracy are understanding the origin of the bad convergence behavior of the pseudo-scalar correlator or computing the values for ratios of the vector correlator on the lattice.  Table 10. α (n f =5) s (m Z ) determination from ratios of charm pseudo-scalar correlator moments R P,n c . The first column shows from which reference the lattice results are taken, the second which ratio has been used, the third corresponds to the central value, fourth to seventh provide the various component of the uncertainty: scale variation (σ pert ), lattice (σ lattice ), charm mass (σ mc ) and gluon condensate (σ np ), which are added in quadrature in the last column (σ total ).

Comparison to previous lattice determinations
In this section we compare the estimates of perturbative uncertainties from the various lattice collaborations, which have a huge impact in the total uncertainty.
Ref. [8] estimates the perturbative uncertainty on α s from the M P,0 c moment by varying the renormalization scales setting 2 GeV ≤ µ α = µ m ≤ 4 GeV. The result 0.1174 (12) is quoted, with a central value in good agreement with our result in table 9, but with a smaller error. If we use a correlated scale variation (taking the charm quark mass as the lowest scale) we obtain 0.1178 ± (0.0011) total , which has a very similar total uncertainty.
Ref. [9] does not perform any scale variation, but estimates the uncertainty by making an educated guess for the unknown O(α 4 s ) term. The 0-th moment is used in the analysis, and value 0.1183(7) is quoted, which agrees well with the corresponding central value in table 9, but has only a third of our uncertainty. Refs. [10] and [11]  , respectively, estimating the truncation error in the same way as in [9]. Our uncertainties turn out to be larger because of our more conservative method to estimate perturbative errors, being our central values a bit larger. By averaging our results which use data from [11], we obtain 0.1162.
Ref. [12] performs uncorrelated double scale variation, taking µ α = µ m ± 1 GeV with 2 GeV ≤ min{µ α , µ m } and max{µ α , µ m } ≤ 4 GeV. The value 0.1177 ± (0.0026) total is obtained from a combined analysis involving M P,n=0,1,2 c and their ratios. The analysis carried out in [12] determines the charm quark mass and the gluon condensate, what may also contribute to a different central value. Our total uncertainty is slightly smaller than theirs, and this could be caused by their different prescription when varying the renormalization scales. If we use an unconstrained grid our total error grows to 0.0029 and 0.0026 for the first two ratios, the latter matching their quoted error.

Final result: fit to charm vector and pseudo-scalar moments
One can go one step further to gain precision by determining α s from a fit to different correlators, which benefits from the fact that the experimental and lattice data are obvi--22 -ously uncorrelated. Since the bulk of the uncertainty comes from truncation errors, it is also crucial to realize that perturbative uncertainties from the vector and pseudo-scalar correlators are completely independent. Hence one can easily construct a χ 2 function in which all uncertainty sources (experimental, perturbative, non-perturbative, and from the charm mass) are combined. The uncertainties due to the gluon condensate and the charm mass are correlated and we implement this into our function. For this analysis we consider two observables: M P,0 c , since results from the lattice collaborations are in excellent agreement, and R V,2 c , which has the smallest error among the vector moment extractions. Adding other moments to the χ 2 would not improve the strong coupling determination, due to the strong correlations among moments from the same correlator. We take an average of all four lattice values, assigning the smallest uncertainty of the four determinations (since perturbative uncertainties overly dominate, this choice has a tiny effect on the fit outcome): M P,0 c = 1.7054 ± 0.0053. On the theoretical side we make a one-dimensional grid in α s , 8 and for each value we scan in the renormalization scales as described in section 4 to determine a central value and perturbative uncertainty for both the vector and pseudo-scalar moments. We make three additional grids: one setting the gluon condensate to zero and other two in which the charm mass is set to 1.3 GeV and 1.26 GeV. These extra grids are used to determine the (correlated) non-perturbative and charm-quark mass uncertainties. Following this procedure one accounts for α s -dependent theory uncertainties in a consistent way.
Our χ 2 function takes the following form where here α s ≡ α The non-diagonal terms in eq. (5.5) arise because of the correlation in the theory uncertainties due to the mass and the gluon condensate, and are written as Upon the minimization of the χ 2 function we find the main outcome of this paper

JHEP03(2020)094
with a p-value of 26%. Neglecting correlations due to the gluon condensate and the charm mass the uncertainty is 0.0013, so once again our treatment is conservative. Changing the value of M P,0 c to a weighted average (with much smaller errors from the lattice) does not change the fit outcome within four digital places. Taking instead the least precise lattice determination lowers the central value by 0.0003, a shift much smaller than the uncertainty.
Our result is compatible with the current word average within 0.65 σ with a smaller central value and a very similar uncertainty. As a sanity check we have also performed a standard weighted average of the individual α s determinations from M P,0 c and R V,2, c , finding α (n f =5) s (m Z ) = 0.1172 ± (0.0014) total , in excellent agreement with eq. (5.8).

Conclusions
In this work we have determined α s from ratios of quarkonium correlator moments. The strategy we follow when taking the ratio is canceling the overall dependence on the heavyquark mass, leaving only a small logarithmic dependence which is damped by two powers of the strong coupling constant. The ratio is re-expanded in α s such that the new series has nice convergence properties. Our determination uses O(α 3 s )-accurate perturbative computations, supplemented with non-perturbative corrections parameterized by the gluon condensate, which is included to O(α s ) precision. We compute the experimental values for the ratios of moments of charmonium, updating the computation performed in ref. [3], and bottomonium, using the results in ref. [5], but, in both cases, keeping α s as a free parameter in the treatment of the perturbative continuum contributions, as well as in the treatment of the light-quark perturbative background in the case charmonium. It is crucial to have full control over the uncertainty correlation among moments, since there are large cancellations of errors when taking the ratios. We also analyze lattice results on ratios of the pseudo-scalar correlator.
We perform a careful study of perturbative uncertainties, and conclude that the renormalization scales of α s and m q should be varied independently in the same ranges as in refs. [3,5]. Furthermore, since the convergence parameter of the series does not show a preferred value, we do not restrict the renormalization scales based on convergence. Instead, we simply impose that the logarithm of ratios of scales is order one. This restriction ensures order-by-order convergence, while keeping uncertainties of moderate size. Our estimates are in general much more conservative than those based on educated guesses for the unknown higher order terms of the perturbative series. In that sense, our reanalyses of lattice data shows that uncertainties are overly dominated by the truncation error, such that more precise lattice simulations will not decrease the total incertitude with the present knowledge on the perturbative series.
Our uncertainty could go further down if the accuracy of the experimental value for the moments of the charm vector correlator increased. On the theory side, since the perturbative error dominates, knowing the moments at O(α 4 s ) would be paramount to obtain smaller uncertainties. We get a value fully compatible with the current world average with a comparable uncertainty. Our results from bottom moments, despite their smaller perturbative errors, are less precise because there is less sensitivity to α s at larger energies even though the experimental values for the ratios are equally accurate. The situation would dramatically improve if the bottom cross section measurements in the continuum reached larger energies, and also with more precise data at threshold. In figure 5, we compare our main result of eq. (5.8) to a selection of previous determinations from a variety of sources. 9 Our analysis can be extended or improved in a number of directions. On the theoretical side, as an additional sanity check on our perturbative error estimate, we could consider alternative expansions, such as taking powers of R X,n q (and re-expanding the resulting 9 See also ref. [76] for a recent extraction based on a strategy similar to the analysis of ref. [77]. -25 -JHEP03(2020)094 series), or using a linearized iterative solution for α s (µ α ) in the same spirit of refs. [3,5]. In the same direction, one could consider analyses in which the ratio of moments is not reexpanded in powers of α s . Finally, since the renormalon associated to the pole mass cancels when taking ratios, an alternative analysis could use directly the pole mass (which appears only in logarithms). In a different direction, QED effects could can be readily accounted for, since final-state photons are included in the cross section measurements. Finally, one could consider fits to α s using all available information (various ratios from charm and bottom correlators) included in a χ 2 function, taking into account all correlations. Similarly, one can perform global fits to α s and the charm and bottom quark masses.
A related analysis could use ratios of bottomonium moments for very large n, using NRQCD predictions for the theoretical moments, which includes resummation of Sommerfeld-enhanced terms, and, in some cases, also the logarithm of the relative velocity of the qq pair, see e.g. refs. [91,92]. In this way, the experimental uncertainty would be reduced to the point of being completely negligible. The question is, of course, by how much the perturbative error grows, and if it would still make the method eligible for a precise determination of the strong coupling. Another approach to exclude the continuum from the experimental moments is using finite-energy sum rules, in which a pinched kernel suppresses quark-hadron duality violations.
In the same spirit of our present analysis, one could consider ratios of masses of bottomonium states. In ref. [38] it was shown that a simultaneous fit for the bottom quark mass and α s was feasible, but the strong coupling was determined with large (perturbative) uncertainties. The origin of the big error was the strong correlation between the two parameters (a similar correlation was found for fits to the Cornell model in ref. [93]). The series resulting from the ratio of two bottomonium bound-state masses (having e.g. distinct principal quantum number n) is free from the pole-mass renormalon and depends only logarithmically on the heavy-quark mass, starting at NNLO. Therefore one a) eliminates the correlation between the two parameters, and b) has a better-behaved perturbative series. Probably one cannot use this idea in charmonium states since masses of mesons with n = 2 are already badly described in perturbation theory.

A Renormalization group evolution for coupling and masses
Our numerical programs use an exact solution to eqs. (2.15), using methods which do not rely on directly solving an ordinary differential equation such as the Runge-Kutta algorithm. For α s , the first line of eq. (2.15) can be integrated as follows: . (A.1) We analytically perform the integration using partial fractions. The root decomposition can always be performed to arbitrary precision once the number of flavors n f has been specified. One is then left with an ordinary equation which cannot be solved analytically beyond LL. Its solution, however, can be easily found numerically e.g. by ordinary methods (as we use in our Mathematica implementation), or using a recursive method with the LL (analytic) solution as seed (used in our python implementation). Solving the MS mass RGE equation is very simple, and as usual one can decouple the α s and m q evolution of eqs. The integration over α s is easily performed with the partial fraction just discussed. In our numerical codes we always use the anomalous dimensions for α s and m q at five-loop order, and the four-loop threshold condition [94][95][96] to relate α We study the convergence properties of the series using V c , the parameter introduced in ref. [5], which is a finite-order version of Cauchy's root convergence test. While in ref. [5] V c was right-away defined on the fit parameter (in that case it was the quark mass), here we apply the method directly to the series. This is justified since the quantity we aim to determine is already the expansion parameter. Writing the series generically as S = i=0 s i with s i ∝ α i s depending on µ α , µ m and m q (µ m ) as ; charm pseudoscalar correlator, 0-th moment (blue), and first two ratios (green and red). For panels (a) and (b) blue, green, and red correspond to the first, second, and third ratios of moments, respectively.