Multiloop contributions to the on-shell-MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{{\mathrm{MS}}}$$\end{document} heavy quark mass relation in QCD and the asymptotic structure of the corresponding series: the updated consideration

The asymptotic structure of the QCD perturbative relation between the on-shell and MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{\mathrm{MS}}}$$\end{document} heavy quark masses is studied. We estimate the five and six-loop contributions to this relation by three different techniques. First, the effective charges motivated approach in two variants is used. Second, the results following from the large-β0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0$$\end{document} approximation are analyzed. Finally, the consequences of applying the asymptotic renormalon-based formula are investigated. We show that all approaches lead to corrections which are qualitatively consistent in order of magnitude. Their sign-alternating character in powers of the number of massless quarks is demonstrated. We emphasize that there is no contradiction in the behavior of the fine structure of the renormalon-based estimates with other approaches if one use the detailed information about the normalization factor included in the renormalon asymptotic formula. The obtained five- and six-loop estimates indicate that in the case of the b-quark the asymptotic character of the studied relation manifests itself above the fourth order of PT, whereas for the t-quark it starts to reveal itself after the seventh order. This allows to conclude that like the running masses, the pole masses of the b and especially t-quark in principle may be used in the phenomenologically-oriented studies.


Introduction
It is well known that the bare masses of quarks in QCD can be expressed through their renormalized finite analogs defined in a particular scheme. In this work we will consider primarily two renormalization schemes, namely the MS-and the on-shell OS-scheme. The latter is used for defining the pole masses of heavy quarks. The relevant renormalization a e-mail: kataev@ms2.inr.ac.ru (corresponding author) b e-mail: viktor_molokoedov@mail.ru prescriptions for masses of these particles have the following form: where m 0,q , M q and m q are the bare, pole and MS-scheme running masses respectively. The renormalization mass constants Z OS m and Z MS m contain the traces of ultraviolet divergences in the form of poles and are represented by the perturbation theory (PT) series in powers of the coupling constant of the strong interaction depending on the scale parameter μ and defined in the corresponding subtraction scheme.
Due to the fact that the masses M q and m q (μ 2 ) are the finite renormalized quantities, their ratio must also be finite. It is convenient to introduce the following relation between the pole and running masses of heavy quarks, also called in the literature as the on-shell-MS mass relation: with the strong coupling constant a s = α s /π defined in the MS-scheme in the Minkowski time-like region. The number of the active flavors n f running inside the fermion loops (the values n f = 4, 5, 6 correspond to the cases of the charm, bottom and top-quarks respectively) is related to the number of the light (massless) n l and heavy (massive) n h flavors by the following way n f = n l + n h . In this work we use the approximation when only one heavy quark is massive i.e. n h = 1 and the rest n l = n f − 1 are massless. The one-loop termt M 1 was calculated a long time ago in [1]. The two-loop correctiont M 2 was analytically computed in [2] and confirmed later in [3,4]. The O(a 3 s ) contribution was evaluated independently by analytical [5] and semi-analytical [6] methods.
In the relation (2) for coefficientst M k (μ 2 /M 2 q ) the transition from the pole mass to the running one can be car-ried out by solving the corresponding renormgroup (RG) equations. After this, it is possible to define the coefficients t M k (μ 2 /m q (μ 2 )). In the normalization point μ 2 In any order of PT the terms t M k can be expanded in powers of the number of n l and n h . Fixing n h = 1 we arrive to the following expansion: In particular, the four-loop coefficient t M 4 is a third degree polynomial in n l : The first two coefficients t M 4,3 and t M 4,2 in (5) were calculated analytically in [7]. Note that the exact numerical expressions of the terms leading in powers of n l in (4) were obtained in [8] up to the ninth order of PT from consideration of the contributions generated by a renormalon-type chain of quark loops inserted into the gluon line which renormalizes the propagator of the heavy massive quark. The last two coefficients in (5) have not yet been computed in the analytical form. However, after the numerical evaluations performed in [9] for the overall term t M 4 at the fixed number n l = 3, 4, 5, the approximate values of the contributions t M 4,1 , t M 4,0 have been obtained with the help of the least squares method (LSM) in [10]. It allows to solve the overdetermined systems of algebraic equations and also to fix the uncertainties of its solutions. The similar expressions for these two terms were also found in [11] by means of a special fitting procedure. It was based on the application of the renormalon calculus of [12][13][14][15][16] and more definitely on the renormalon asymptotic formula for coefficients of the relation between the pole and running masses of heavy quarks originally derived in [15,16].
Later on the evaluation of the t M 4 -coefficient was done in [17] with higher precision than in [9] and for a much larger number of flavors in the range 0 ≤ n l ≤ 20. The central values of terms t M 4,0 and t M 4,1 , being also extracted in [17] from the fitting of the numerical results for t M 4 at the fixed number of n l , agree with the ones obtained with the help of the LSM and presented in the "Note added" of [10] and in [18].
Despite the apparent smallness of the four-loop corrections to the on-shell-MS heavy quark mass relation, its knowledge is important from both phenomenological and theoretical points of view. Indeed, the asymptotic nature of the PT series 1 for the relation between the pole and running masses, which is governed by the dominant u = 1/2 infrared renormalon contributions to the Borel image of this relation (discovered in [13,14]), leads to the factorial growth of the coefficients t M k at the large orders k. This fast increase is associated with the sensitivity of the pole mass to small momenta, due to which it suffers from the large perturbative corrections [14,15]. On the contrary, the running mass, defined within the MS-scheme, depends on the ultraviolet (UV) subtraction of divergences only and, therefore, does not contain the infrared (IR) renormalon contributions. In this regard, it is very important to know in the specific physical studies when the asymptotic behavior will begin to manifest itself in the definite cases of the charm, bottom and top quarks.
As follows from the results of [2][3][4][5][6] for the c-quark the asymptotic behavior of the on-shell-MS relation reveals itself in the rather low orders of PT, namely in the second or third order (depending on the normalization point). Therefore, in the modern high-precision studies it is more preferable to use the concept of its running mass with the value extracted e.g. in [20][21][22][23][24][25].
In the case of the b-quark the first traces of the asymptotic structure of the ratio s ) level [9,17]. However, for an unambiguous response to the question about the number of order of PT starting from which the asymptotic behavior will manifest itself it is necessary to know the value of the correction of the fifth order. But already from the available data, it follows that unlike the running mass at the four-level the pole mass of the bottom quark should be used with care. The values of m b (m 2 b ) were obtained at the N 3 LO level as the final results of the QCD analysis of the properties of Υ system within the static potential studies (see e.g. [23,[25][26][27][28]), the QCD sum rules (see e.g. [20,21,29]) and of the production cross-section of the bb-quarks in the e + e − collisions [30]. Also worth mentioning the recent results of the lattice QCD determinations [31]. These lattice results are stimulating a more careful study of the existing uncertainties in the four-loop on-shell-MS mass relation for the b-quark.
For the t-quark the situation is even more intriguing. The definite results of the experimental analysis of Tevatron and LHC data are expressed through a Monte-Carlo t-quark mass, which may be related (though with process-dependent uncertainties) to its pole mass (for the detailed consideration see e.g. [32][33][34]). The average PDG (20) value of this important quantity, obtained from the recent LHC measurements and the updated Tevatron analysis, is M t = 172.76 ± 0.30 GeV [35]. For comparison, recently the LHC value of the pole mass of the t-quark with the thorough estimates of the various types of uncertainties was obtained in [36] and reads M t = 171.1 ± 0.4(stat) ±0.9(syst) +0.7 −0.3 (theor) GeV. Note that in the process of getting the top-quark mass values the question about inaccuracies of different Monte-Carlo programs used for analyzing Tevatron and LHC data has become more vivid. This problem is still under careful examinations (see [37,38] and the reviews of [33,34]). The arising inaccuracies should be compared with other theoretical errors, which enter into the determinations of both running and pole top-quark masses [24,39,40]. Moreover, the uncertainties, contributed by at least the first not yet computed highorder correction to the relation between the pole and running masses, are also of interest [38,41]. The study of these effects will be continued in this paper using several approaches for estimating high-order QCD corrections to the on-shell-MS mass relation.
Our main aim is to analyze the asymptotic structure of the perturbative series for the ratio M q /m q (m 2 q ) at the O(a 6 s ) level. To get a feeling for what may be the values of the five-and six-loop corrections to this ratio, we estimate them using three distinct techniques. After this, we restore the general n l -dependence of these estimates (the previous definite results on this topic are presented in brief in [42][43][44]) and demonstrate its sign-alternating character in n l .
The outline of our studies is as follows. In Sect. 2 we present the current known four-loop corrections to the onshell-MS mass relation for the particular case of the SU (3) color gauge group. Here we especially emphasize the appearance of the contributions proportional to powers of π 2 -terms to the analytical expressions for coefficients of the ratio (3) starting to manifest itself from the two-loop level and originating from calculation of Z OS m in the Minkowskian on-shell subtraction scheme (the first emergence of a π 4 -term in Z MS m occurs at the O(a 4 s ) level only). In Sect. 3 we use the Källen-Lehmann type dispersion relation for the "effective" spectral function, defined in the Euclidean domain for energies, to model the on-shell π 2terms contributing to Z OS m by the analytical continuation π 2 -effects arising upon the transition from the Euclidean to Minkowskian region.
In Sect. 4 we apply the approach proposed in [45] and extended in [46] to estimate five-and six-loop corrections t M 5 and t M 6 with partial incorporation of the π 2n -contributions being mentioned above. This approximate procedure is based on the effective-charges (ECH) method [47] and on the concept of scheme-invariants [48].
Section 5 is devoted to the study of the consequences following from the results of [8], where the exact numerical values of the contributions leading in powers of n l to the coefficients t M k were computed from consideration of the diagrams containing an insert of a chain of quark loops into the single gluon line, renormalized the massive quark propagator. Note that within the Naive-Nonabelianization (NNA) procedure utilized by us, this leads to the sign-alternating n l -structure of the five-and six-loop PT corrections. Section 6 is dedicated to the investigation of the O(a 5 s ) and O(a 6 s )-estimates found with help of the asymptotic renormalon-based formula for coefficients of the on-shell-MS relation, which was previously studied in [15,16,41,[49][50][51]. Herewith, we consider two variants for fixation of the normalization factor included in this factorial formula (for details see [41] and [50]). We demonstrate that using both these ways one can obtain the sign-alternating structure of the five-and six-loop coefficients in (3). This fact is in full agreement with the outcomes following from the application of the ECHmotivated method and the large-β 0 analysis. Here we especially emphasize that in contrast to the results of our previous works on this topic [42][43][44] the sign-alternating structure of the renormalon-based estimates is observed upon attraction of more detailed information on the normalization factor of the renormalon asymptotic formula.
In Sect. 7 we briefly summarize all the main our results presented in the previous sections and consider the numerical impact of the estimated O(a 5 s ) and O(a 6 s ) terms on the behavior of the on-shell-MS relation for real heavy quarks. We show that the application of all methods employed by us leads to the results which are consistent with each other in order of magnitude (on average with a factor two).
For clarity in Appendix A of this paper we set out the key points of the LSM, define the way of finding the LSMsolutions for the terms t M 4,0 and t M 4,1 and their uncertainties. Note here that as follows from the studies of [10,18] these solutions of the overdetermined system of algebraic equations are stable under a change not only in the number of n l -equations being considered, but also in the number of unknowns involving in this system. In order to consider the possible differences in the structure of the perturbative series in QCD and QED in Appendix B we compare the behavior of the PT series for the relation between the pole and running masses of the heavy quarks in QCD with the corresponding one for the charged leptons in QED at the four-loop level.

The on-shell-MS heavy quark mass relation: available analytical perturbative QCD results
Consider first the relation (3) between the pole and running heavy quark masses normalized at the scale μ 2 = m 2 q . It is known that the heavy quark pole mass M q is defined in the on-shell scheme as a pole of the renormalized heavy quark propagator in the Minkowski region. In turn, the scale evolution of the MS-scheme heavy quark running mass is first defined in the Euclidean domain since the calculations of the corresponding master integrals for Z MS m are also performed in the Euclidean region: This relation may be transformed to the Minkowski region by replacement Q 2 → s. After this, it is possible to fix the Minkowskian scale μ 2 = m 2 q and to define m q (m 2 q ). The dependence of the QCD expansion parameter a s (μ 2 ) and of the running quark mass m q (μ 2 ) on the renormalization scale μ 2 is determined by the following RG equations: where β(a s ) and γ m (a s ) are the QCD β-function and the anomalous mass dimension. In our further consideration we use their MS-like scheme expressions. The one-and two-loop coefficients β 0 and β 1 of the QCD β-function were computed analytically in [52][53][54][55][56] respectively. The symbolical expressions of the scheme-dependent three-and four-loop coefficients β 2 and β 3 are known from calculations performed in [57][58][59][60] correspondingly. The coefficient β 4 was obtained in analytical form in the SU (3)-group [61] and confirmed in [62,63] by computing this term in the general SU (N c ) gauge group. Note that in the process of these calculations the Euclidean contribution, proportional to the ζ 4 = π 4 /90 Riemann function, is appearing for the first time. For our purposes it is convenient to present these coefficients β n in terms of the number of massless flavors n l = n f − 1. In the case of SU (3) color gauge group their numerical expressions have the following form: The first scheme-independent coefficient γ 0 of the QCD anomalous mass dimension function of Eq. (8) was presented in [1]. Its two-, three-and four-loop expressions were analytically computed in [1,[64][65][66][67][68] correspondingly. The coefficient γ 4 of the fifth order was evaluated in case of the SU (3) color gauge group in [69]. This analytical result had been confirmed later on in [70] upon more general calculations performed in the SU (N c )-group. It should be stressed that the Euclidean contributions ζ 4 being proportional to π 4 are arising in the QCD expression for γ m beginning from the four-loop level (see [67,68]), whereas the functions ζ 6 proportional to π 6 are starting to reveal themselves at the fiveloop level. 2 The numerical values of these coefficients are: Note that all renormalized quantities, which enter into ratio (3), are self-consistently defined in the Minkowski region of energies. In particular case of the SU (3) group the analytical contributions to the first four coefficients of Eq. (3), expanded in powers of n l (see Eq. (4)), follow from the calculations of [1,2,5,7] and read  3 where Li n (x) = ∞ k=1 x k k −n is the polylogarithmic function.
As the result the two-and three-loop coefficients of the ratio M q /m q (m 2 q ) have the following numerical form: The numerical n l -dependent expression for the O(a 4 s ) term t M 4 is known at present with high enough accuracy. We combine here the results of the analytical (11g-11h) [7] and semi-analytical computations [17] with the LSM-solutions [18] for the constant t M 4,0 and linearly dependent on n l term t M 4,1 with their LSM-uncertainties. This leads to the following expression (see Appendix A): Note that for our purposes to study the asymptotic structure of the on-shell-MS mass relation the uncertainties included in (12c) are not important and we can neglect them.
Unlike the coefficients of the QCD β-function and the anomalous mass dimension the results (12a-12c) clearly demonstrate the sign-alternating pattern in n l . It is interesting to note that this computational fact is consistent with the theoretical renormalon-inspired large β 0 -expansion [8,12].
We now return to the discussion concerning the analytical structure of certain contributions to the formulas (11b-11h). It is worth emphasizing that the second, third and fourth coefficients t M 2 , t M 3 , t M 4 contain the π 2 -terms typical to the Minkowskian on-shell subtraction scheme, while the additional π 4 -contributions are emergering in the results of threeloop calculations and beyond. We expect the appearance of the π 6 -terms in the structure of analytical expressions for the yet unknown coefficients t M 4,0 and t M 4,1 . 3 Comparing the analytical structure of the O(a 4 s ) perturbative QCD corrections to the ratio M q /m q (m 2 q ) and to the QCD anomalous mass dimension γ m (a s ), dictated by the pattern of the quark mass renormalization constant Z MS m , one can conclude that only the π 4 -contributions, entering into the t M 4,2 and t M 4,1 , may contain the admixture of the typical Euclidean ζ 4 -terms, first appearing in the four-loop contributions to γ m (a s ), which are proportional to n 2 l and n l [67,68]. Other contributions to the coefficients t M 2 − t M 4 , proportional to powers of π 2 , arise 3 This expectation is already supported by the recent QED analytical results of t M 4,0 in [72].
from computations of high-order corrections to the renormalization constant Z OS m defined in Eq. (1) in the Minkowskian on-shell subtraction scheme. In next section we try to build an analogy between these typical on-shell scheme π 2 -contributions and the "kinematic" effects proportional to powers of π 2 in the perturbative QCD expressions for the Minkowskian physical quantities, initially evaluated in the MS-scheme in the Euclidean domain. The substantial role of these effects has been demonstrated in the number of works on the subject (see e.g. [73][74][75][76][77][78][79][80]).
3 Is it possible to link the on-shell π 2 -contributions and "kinematic" π 2 -effects?
To understand whether it is possible to draw the analogy between contributions proportional to powers of π 2 to the ratio M q /m q (m 2 q ), which are defined in the Minkowskian region and the "kinematic" π 2 -terms, arising in the PT coefficients of the Minkowskian RG controllable physical quantities in the MS-scheme and associated with the analytical continuation effects from the Euclidean to Minkowskian domain, we will follow the path treaded in [46] and used later on in [81]. For this goal, we consider the Källen-Lehmann type dispersion representation, 4 which allows to simulate the appearance of these "kinematic" terms: Here the model spectral function T (s) is determined in the Minkowski region 5 as: In this perturbative expression m q (s) is the MS-scheme running mass of heavy quark, normalized at the scale μ 2 = s in the time-like region and t k are the dimensionless coefficients of this spectral function. 6 One of the basic ideas of the work [46] consists in a fact that at s = m 2 q the coefficients t k in Eq. (14) are assumed to be equal to the corresponding onshell scheme coefficients t M k of the heavy quark mass relation Substituting the expression (14) into relation (13) one can arrive to the perturbative representation for the Euclidean function F(Q 2 ) with coefficients f E k related to t M k by the following way where contributions Δ k are the "kinematic" terms, which reflect the analytic continuation effects. Note that from the point of view of the first principles of the theory of dispersion representations the model equation (13) is not completely substantiated. Indeed, within PT it should contain the subtraction constant, which is related to the theoretical ambiguities in the low-energy region, discussed in [80,82] upon a study of the dispersion representations of the Green's functions for the scalar quark and gluon currents. In this regard, it would be more consistent to consider the model subtracted dispersion relation written down for the function F(Q 2 ) − F(0). However, below we will show that the perturbative estimates for coefficients of the ratio M q /m q (m 2 q ), obtained with help of the expression (13), yield the quite reasonable predictions of the asymptotic behavior of this ratio and agree with applications of the renormalon-motivated calculus (with a factor of order 2).
Keeping in mind the aforesaid discussions, Eqs. (13)(14), remark on equality of coefficients t k and t M k at s = m 2 q and taking into account the inverse integral representation for the function T (s) 7 one can obtain the following approximate representation for the pole and MS-scheme masses of heavy quarks [46]: Using now Eqs. (13)(14)(15) we can fix the explicit form of the "kinematic" contributions Δ k in (16) up to the sixth order of PT. Far enough from the regions of manifestation of the heavy quark threshold effects the differential system of RGequations (7)(8) in the time-like region can be rewritten in 7 By analogy with the dispersion relation between the e + e − annihilation R(s)-ratio and the Adler D(Q 2 )-function, the integration contour on the plane of complex variable z lies in the region of analyticity of the integrand (here function F(Q 2 ) is the analog of the Adler function). the following integral form in the O(a 6 s ) approximation: Substituting solutions of this system into function T (s) in Eqs. (13)(14) we get the following integrals which are equal to: where l = ln(μ 2 /s) and L = ln(μ 2 /Q 2 ). Fixing further μ 2 = Q 2 we find the explicit expressions for the terms Δ k .
In the recurrent form they read: The terms Δ 1 − Δ 4 agree with the ones, obtained previously in [46]. The expressions for Δ 5 and Δ 6 are new. One can see that the six-loop contribution Δ 6 does not contain yet unknown coefficients β 5 and γ 5 . They are included only in terms depending linearly on ln(μ 2 /s), which due to Eq. (19) vanish automatically in the Euclidean renormalization point Taking now into account the relation (16) and numerical expressions for the coefficients of β(a s ), γ m (a s ) and t M k , given in Eqs. (9a-10d), (12a-12c), we arrive to the following numerical n l -dependent results for the terms Δ k : where in the expression for Δ 6 -contribution we have neglected the relatively small mean square errors following from computations of the coefficient t M 4 [17]. Worth emphasizing that despite the non-regular sign polynomial structure of the coefficients of the QCD RG functions β(a s ) and γ m (a s ) (9a-10d), the analogous expressions for contributions Δ k respect the alternation of signs in powers of n l that is typical to the two, three and four-loop coefficients t M k .
Their numerical values for the specific number of massless flavors are presented in Table 1.
The outcomes of Table 1 demonstrate the significant growth of terms Δ k with increasing of an order k of PT. This effect is determined by two factors. The first of them is related to the factorial rise of the coefficients t M k included in the definition of the contributions Δ k (20b-20e) (see Sect. 6 of this paper, where the renormalon-based asymptotic formula for t M k is discussed). The second one is partially associated with the considerable factorial growth of the constant terms appearing in r.h.s of Eq. (19) upon the integration of the RG-logarithms with various degrees. Indeed, as was shown in [83] the dimensionless analog of integral (19) with arbitrary degree n has the closed form where the variable x = s/Q 2 and terms ζ 2 = π 2 /6, ζ 4 = π 4 /90, ζ 6 = π 6 /945 may be explicitly restored in r.h.s of Eq. (19). Since 1 < ζ 2 p ≤ ζ 2 < 2 for any p ∈ N, then at even n the integral (22) is factorially growing. As a result, the constant terms in r.h.s of Eq. (19), which enter in the contributions Δ k , are factorially increasing with order of PT as well. Moreover, matching Eqs. (20a-20e) with (22) we conclude that the contribution to the even-order term Δ 2 p leading in powers of π 2 behaves itself by the following way for any p ∈ N: Thus, we conclude that the overall Minkowskian "kinematic" π 2 -effects are indeed not small. Moreover, the values of Δ 2 , Δ 3 , Δ 4 are comparable with the corresponding coefficients t M k (n l ) (see Eqs. (B.13a-B.13c) in Appendix B). In this regard and in view of our assumption that these fast growing "kinematic" effects may model the π 2 -contributions to the high-order coefficients of M q /m q (m 2 q )-ratio, typical to the on-shell renormalization scheme, we note that it is really worth to treat these special terms with care.

The effective charges-inspired estimates
Let us first study the application of a variant of an RGinspired approach for estimating high-order perturbative corrections to a physical quantities being formulated and developed in [45]. This approach is based on the effectivecharges (ECH) method [47]. In the work [46] it was first adapted to the quantity F(Q 2 )/m q (Q 2 ) defined in the Euclidean region. Since here we consider the case of n l massless flavors, running inside the fermion loop insertions of a self-energy operator renormalizing the massive heavy quark propagator, then the coefficients f E n in Eq. (15) do not depend on masses. In this approximation the perturbative expression for F(Q 2 )/m q (Q 2 ) is also independent on mass parameter. Therefore one can use directly the methods described in [45,47]. The corresponding ECH coupling constant a e f f s (Q 2 ) may be defined as: where The coefficients of the ECH β-function for a e f f s (Q 2 ) are expressed through scheme-independent combinations [48] of the higher order PT contributions φ k (25) and β k of the MS-scheme β-function. At the four-loop level these combinations have already been applied for determination of the ECH β-function of the static potential in the QCD [84]. Here we present the explicit expressions for six coefficients of the corresponding ECH β-function, which is governing the : Our further analysis is based on the theoretical studies described in [45,46,81]. Their essence was as follows: if one put β e f f 2 ≈ β 2 , then from Eq. (26b) one can get that afterwards. Estimates of this type were made in [46,81] to fix the numerical value of the term t M 4 for the cases of the charm, bottom and top-quarks which was still unknown at that time.
An admissibility of the approximation β ≈ β k , we will get the estimates of t M k -terms directly without additional evaluation of the "kinematic" Δ k -corrections (t M, EC H direct k stands for these estimates). Nevertheless, these terms will include the π 2 -contributions typical to the on-shell scheme. Indeed, in the estimates t M, EC H direct k of k-th order these π 2 effects are contained in the known analytical Minkowskian coefficients of (k − 1)-th order and lower. This approach will be considered in more details below.
The  Table 2.
One can see from data of Table 2 that both variants of the ECH-motivated method give quite good approximations for the three-and four-loop coefficients of the ratio M q /m q (m 2 q ) 8 (apart from the non-physical case of n l = 8 for the ECH approach, where the estimation differs from the genuine value by a factor over 3). Indeed, both these implementations predict not only the correct signs for the coefficients of the O(a 3 s ) and O(a 4 s ) orders but also yield the estimates whose values are rather close to the expressions having been calculated exactly.
Let us now probe the n l -dependence of the estimated coefficient t M 3 with three unknowns t M 3,0 , t M 3,1 , t M 3,2 with help of three physical data points 3 ≤ n l ≤ 5. Solving the corresponding system of equations we gain the following expansions: The sign-alternating structure of the n l -expanded expression (28b), gotten in the Minkowskian region directly, is consistent with result (12c) of the explicit diagram-by-diagram calculations. However, the ECH approach, applied in the Euclidean region, leads to the different sign of the leading cubic term in (28a). To clarify this observation we note that its value almost coincides modulo with n 3 l -contribution to Δ 4 in Eq. (21c). This means that the cubic coefficient in f E 4 = t M 4 + Δ 4 is close to zero. For a more detailed study of this fact we present Fig. 1, where the obtained expansions It is seen from Fig. 1 that the relative uncertainties of the ECH-direct approach are stable to the changes of n l . However, this is not true for the Euclidean ECH method applied for estimation of the term t M 4 . Indeed, in this case the relative error varies in a wide range from 25% at n l = 3 to 5% at n l = 6. Moreover, at n l = 8 we observe a mismatch with the exact result by a factor over 3 (see Table 2). Therefore, we conclude that at the relative errors of about 25%, one should not trust the almost zero estimated value of the cubic n l -dependent term in f E 4 coefficient. Thus, the errors of this order are quite satisfactory while getting the estimates of the term t M 4 (or f E 4 ) at the fixed number of massless quarks but they turn out to be unsatisfactory for the study of more subtle effects of its flavor dependence. Therefore, we infer that the mismatch of the sign of n 3 l -term in Eq. (28a) to the true one is a rather accidental fact related to the instability of the uncertainties being discussed above. In view of this we do not consider the positive sign of this cubic coefficient as a violation of the indication of the sign-alternating structure of the n l -expanded ECH-based estimates of t M 4 -term following from (28b).
The acceptable agreement of the ECH-estimates of the coefficients t M 3 and t M 4 with the results of explicit calculations at the fixed number of n l allows us to regard both variants of the ECH-inspired method as satisfactory estimating procedures. Therefore, we will apply these two realizations to evaluate the unknown contributions of the fifth and sixth orders of PT to the on-shell-MS heavy quark mass relation as well.
Our further studies of the ECH method, initially applied to the Euclidean physical quantities, will contain the following steps: 1. At the first stage, using the explicit expressions for the known terms t M 1 − t M 4 (11a), (12a-12c) and adding to them π 2 -effects of the analytical continuation (21a-21c), we find the Euclidean ing β e f f 5 ≈ β 5 in (26e) and using the numerical expression for f E 5 , obtained at the previous stage, 9 we primarily estimate the value of f E 6 -contribution and then, taking into account (16) and (21e), evaluate the value of t M 6correction to the ratio M q /m q (m 2 q ).
In stages 2 and 4, we get the following five-and six-loop coefficients f E 5 and f E 6 of the Euclidean quantity (15): Table 3 The estimates of the coefficients t M 5 and t M 6 , obtained within two variants of the ECH-approach (ECH and ECH direct), the large β 0expansion (FL and FL, M → m) and the asymptotic IRR-based formula (r-n) with the normalization factor N m taken from the direct analysis and consistent with the renormalon sum rule approach (for details see 6.1 and 6.2). The fourth and fifth columns correspond to various choices of the initial scales in the results of [8]  +5 Utilizing these expressions and Eqs. (21d-21e) for Δ 5 , Δ 6 that model the "kinematic" π 2 -terms, we estimate the ECH-values of the coefficients t M 5 , t M 6 at the fixed number of massless flavors given in Table 3. There we also present the estimates of the same coefficients, obtained by us from the application of the large-β 0 expansion to the renormalonchain contributions [8] to coefficients of the on-shell-MS mass relation (see 4-and 5-th columns of Table 3) and with the help of the infrared renormalon (IRR) asymptotic formula [16,41,49], used recently in [25,50,85] (see 6-th column). The details of these our analyzes will be discussed below.
The data from Table 3 demonstrates that at the physical values of n l the obtained estimates agree with each other at the level of factor two. However, the theoretical uncertainties increase drastically starting from the non-physical sector n l ≥ 6. Indeed, on the contrary to the results following from the application of the NNA procedure to the outcomes of [8] (see columns 4 and 5) and from the renormalon asymptotic formula [16,41,50] (see column 6), the ECH-based estimates at n l = 7, 8 take the negative values. The similar signchanging feature also reveals itself in the renormalon studies from n l ≥ 9 (see e.g. the analysis of [27] and [41]). Moreover, the indication of the sharp growth of the uncertainties in the non-physical sector of n l also follows directly from the renormalon studies [50] (see discussions in Sect. 6 below). Therefore, in this work we restrict ourselves by the consideration of the values of n l from the range 3 ≤ n l ≤ 8. This number of data points is definitely enough to investigate the flavor dependence of the six-loop coefficient t M 6 . Let us study the n l -dependence of the ECH coefficients whose numerical values at the fixed number of n l are given in Table 3. As follows from (4) the five-loop contribution t M 5 is the fourth degree polynomial in n l , namely t M 5 = t M 5,4 n 4 l + t M 5,3 n 3 l + t M 5,2 n 2 l + t M 5,1 n l + t M 5,0 . It contains five unknown terms t M 5,4 − t M 5,0 . Therefore, in order to get their numerical values we will use five equations only which follow from the data of Table 3 The numerical solution of (30) leads to the following expression: In the case of the coefficient of the sixth order of PT the similar consideration at 3 ≤ n l ≤ 8 yields: Both expansions (31a-31b) have the sign-alternating structure in powers of n l , which is supported by results of the large-β 0 analysis [8]. Thus, the ECH-motivated method, applied initially in the Euclidean domain and supplemented by the analytical continuation π 2 -effects, leads to the n ldependent structure of the terms t M 5 and t M 6 , which is similar to the ones observed for the exactly calculated corrections t M 2 − t M 4 given in (12a-12c). Let us consider what will happen with the expressions (31a-31b) if one fix in them t M 5,4 ≈ 0.9 and t M 6,5 ≈ −1.5, following from the exact numerical computations of Ref. [8]. These results were obtained there from the consideration of the subset of the specific renormalon-chain diagrams for the heavy quark propagator. Since in this case one less n l -dependent term should be defined, we exclude from the analysis the data points n l = 7(8) upon the estimation of the 10 The square matrix in l.h.s. of (30) is the Vandermonde matrix. It possesses the interesting mathematical properties: the elements of its each row are the terms of a geometric progression and its determinant is equal to Δ = 0≤i< j≤k ((n l + j) − (n l + i)) = 0≤i< j≤k ( j − i). Here the number of massless quarks varies from n l to (n l + k), k ∈ N. flavor dependencies of the coefficients t M, EC H 5 (6) correspondingly. This leads to the insignificant changes in all coefficients of the expressions (31a-31b) with keeping their signalternating character. 11 Using numbers shown in the third column of Table 3, we also obtain the approximate n l -dependence of the O(a 5 s ) and O(a 6 s ) coefficients of the on-shell-MS heavy quark mass relation within the ECH-motivated approach, applied directly in the Minkowskian region:  Table 3), both realizations of the ECH method predict not only the sign-alternating structure of these corrections in powers of n l but also lead to the values of the separate n l -dependent terms close in magnitude. Note also that the fixation of the known terms leading in n l does not substantially affect the values of other coefficients in (31c-31d). 12 In the next sections we will compare these results with the similar ones which follow from the large-β 0 approximation [8] and from the asymptotic renormalon formula [15,16] subsequently improved in [41,49,50].

The consequences of the leading renormalon chain calculations
Before the analytical computations [5,7] of the leading O(a 3 s ) and O(a 4 s ) n l -contributions to the coefficients t M 3 , t M 4 (see (11d) and (11g)), these terms were evaluated numerically in [8]. These results follow from calculations of the leading renormalon-type contributions, generated by a chain of the fermion loop (FL) insertions into the gluon line, renormalizing massive heavy quark propagator. The outcomes of Ref. within the large-β 0 expansion and get their flavor dependencies. Since the terms leading in n l do not depend on μ 2 , we will consider the five and six-loop estimates in two scale normalizations, namely μ 2 = m 2 q and μ 2 = M 2 q with its subsequent transition to the running mass.
Using the results of work [8] and assuming the normalization at μ 2 = m 2 q , we get the following expansions: Next, presuming that the initial normalization point is fixed on the pole mass and then it is shifted to the running one, we find the following analogs of (32a-32b): These expressions demonstrate that the FL-method supplemented by the NNA procedure gives the strict alternation of signs in the polynomial flavor decomposition of the terms t M 5 and t M 6 . This feature is the direct consequence of the application of the large β 0 -expansion. Indeed, in this approximation the k-th term t M k is proportional to β k−1 0 (n l )-factor, where the first coefficient of the QCD β-function is defined in (9a). Therefore, this approach will always lead to the signalternating n l -structure of the estimated corrections in all orders of PT. Moreover, this statement does not depend on the normalization point. Thus, the FL-approach supports the results of the ECH-method presented by us above.
Note that the specific n l -dependent terms in (32c-32d) are smaller than the corresponding ones in (32a-32b). Herewith, the latter are closer to those found with help of the ECHmotivated method in (31a-31b), (31c-31d). The numerical values of the corresponding FL-estimates at the fixed number of n l are presented in the fourth and fifth column of Table 3.

Renormalon-based estimating procedure
Let us now move on to the consideration of another approach for estimation of the high-order corrections to the relation between the pole and MS running masses of heavy quarks based on the renormalon analysis. It is known that the ratio M q /m 2 q (m 2 q ) contains the linear infrared renormalon (IRR) contributions, which lead to the rather strong factorial increase of the coefficients in this asymptotic PT series [13,14,16]. This fast growth of t M k -terms is governed by the leading u = 1/2 IRR pole in the Borel image of the ratio being discussed. The study of the behavior of t M k -coefficients in the renormalon language results in the following asymptotic formula derived in [15,16,41,86]: where Γ (x) is the Euler Gamma-function, b = β 1 /(2β 2 0 ) and the values of the sub-leading coefficients s k , evaluated in [16,41,49,50], are presented below. In the finite order of PT the factor N m depends on n l and k. Note that our normalizations and notations for the coefficients of the QCD β-function (9a-9e) differ from those used in [14,15,41,49] upon studying the formula (33). Indeed, in these works the analytical expression for the first coefficient of the RG β-function of the SU (3) QCD is defined as β 0 = 11/4 − n l /6, while we are using β 0 = 11/4 − (n l + 1)/6 (see (9a)). To coordinate these notations and use directly the asymptotic formula (33) we need to perform a shift n l → (n l − 1) in Eqs. (9a-9e). In this section we will work in these designations.
The corresponding expressions for the coefficients s k read: One should also mention that another (recurrent) way of obtaining the formula (33) was considered in [85]. It was based on the fact that the leading renormalon contribution to the relation between the pole and running masses of heavy quarks is independent on the MS-scheme mass m q (for details see [15] and [50,87]). On the one hand, this fact allows to use the approximate relation d M q /dm q (m 2 q ) ≈ 1. On the other hand, the RG-based form of this derivative can be obtained from Eq. (3) with taking into account the running of the coupling constant. Matching the unit to this RG-based expression one can obtain the recurrence relation which results in to the factorial formula (33).
The ways to fix values of the normalization factor N m in a specific finite order of PT are different in various works. For instance, in [27,41] they were extracted from juxtaposition the results of the diagram-by-diagram calculations of the coefficients t M k with the ones being rewritten in the asymptotic form (33). In the works [27,49,50,88,89] the approximate analytical expressions for N m were obtained out of the analysis of the behavior of the Borel image of the ratio (3). In this paper we will utilize both these ways for fixing of N m -factor.
The extraction of its accurate values for a certain number of n l is extremely important for the investigation of the subtle effects of the n l -structure of the coefficients t M k . As we will show further the change in the accuracy of the used values of N m from two accounted decimal digits to three ones affects considerably this flavor structure. This understanding is supported by the existence of the renormalon sum rule expressions for the normalization factor N m [49,89]. 13

The direct analysis
In the first approach, which we will call the direct one, the numerical O(a 4 s ) values of N m were found in [41] from comparison the results of the explicit four-loop computations [17] with the numerical expressions that follow from the large k renormalon-based expectations (33). The results of this analysis are presented in the second row of Table 4 and are labeled as direct. Note that they are consistent with the ones obtained in [27]. In the third line of this table the numerical values of N m , defined within the renormalon sum rule approach [50] (see discussions below), are given.
The comparison of the direct four-loop values of N m , presented in Table 4, with the three-loop ones, extracted in a similar way in [41] using the O(a 3 s ) results of calculations [5,6], demonstrates that at least at the fixed physical number of flavors the values of N m -factor are rather stable to the transition from one order of PT to another. Indeed, at n l = 3, 4, 5 the relative difference between them does not exceed 15%, but increases substantially for the non-physical flavors. This observation yields us grounds to expect that at least at 3 ≤ n l ≤ 5 the application of the renormalon-based asymptotic formula (33) with N m , taken in the four-loop approximation, will lead to the quite acceptable estimates for terms t M 5 and t M 6 . However, to investigate the n l -structure of these terms in our further analysis we will also consider 6.2 The renormalon sum rule approach Consider now the values of N m , given in the third line of Table 4 and fixed in work [50] by the renormalon sum rule: where the expressions for the terms S k contain the high-order coefficients of the QCD β-function and the corrections t M k of the ratio (3). This formula of [50] was obtained from the consideration of the Borel image of the relation between the pole quark mass and its MSR-mass: The concept of the low-scale MSR-scheme quark mass m MSR q (R) is related to the MS-mass but may be evolved to the renormalization scale R below mass of the heavy quark being considered (for details see [50,87]). The MSR-mass coincides with the MS running mass at R = m q (m q ). An important feature of (37) is that unlike the relation between the pole and running quark masses this expression contains the linear term in R. This linear dependence allows to get the Borel transformation function of (37) in analytical form and the analytical expression (36) for normalization factor N m in particular [50]. Moreover, as was mentioned in this cited paper, this feature permits to speed up the convergence rate of the series (36) compared to the analogous one, obtained in [49,90] as a residue of the Borel transformation function of the ratio M q /m q (m 2 q ) in the leading IRR pole u = 1/2. Note that the study of the effects of other IR and UV renormalons was considered recently for e.g. in [51].
The corresponding uncertainties of N m -factor presented in Table 4 were extracted in [50] from the scale variation of R, which model an effect of the missing higher orders of PT to the ratio (3). These uncertainties increase with the growth of n l . Indeed, at n l = 3 the relative error of N m makes up 2.3% and at n l = 7 it is about 23.2%. At n l = 8 the absolute error rises sharply and already exceeds the central value of N m . This is one more argument why we do not include in our analysis the non-physical sector with n l ≥ 8 upon defining the flavor dependencies of the terms t M 5 (t M 6 ). It is apparent from the data of Table 4 that within the error bands the values of N m , obtained in the renormalon sum rule approach, are consistent with those found by the direct way. Therefore, the expansions (35a-35b) may be valid for this method as well. 15 At the fixed number n l = 3, 4, 5 the numerical estimates t M, r −n 5 and t M, r −n 6 , given in Table 3, are in good agreement with the ones obtained in [50]. Moreover, in the case of the b-quark our five-loop estimates are consistent with those obtained in [25] from the global fits of the energy of the QQ bound states.

Discussion
In this section we briefly summarize all main results obtained above and consequences following from them. Accordingly to the outcomes presented in the previous sections, there are indications that the flavor structure of the five-and six-loop coefficients in the relation between the pole and running masses of heavy quarks has the sign-alterna-ting character in powers of n l , analogous to the one observed for the two-, three-and four-loop coefficients having been calculated exactly. It should be noted in passing that this behavior agrees with the results of the large-β 0 expansion.
Comparing the magnitudes of the estimated leading terms in n l , obtained by us in two variants of the ECH-method and with help of the IRR-based formula, we conclude that the results (31c-31d) in more extent corresponds to the known coefficients t M 5,4 and t M 6,5 (see also the related expressions presented in footnote 12).
For more clarity, we accompany the estimated flavor dependencies (31a-31b), (31c-31d), (32a-32b), (32c-32d), (35a-35b) of the coefficients t M 5 and t M 6 with the corresponding plots, presented in Fig. 2. Figure 2 shows that the estimates, obtained in this work with help of all different methods considered by us, are qualitatively consistent with each other (on average with a factor two) and lead to the rather similar flavor structures of the coefficients t M 5 and t M 6 . Let us now consider an impact of the O(a 5 s ) and O(a 6 s ) QCD estimations gotten within all studied approaches on the behavior of the PT series for real heavy quarks in more detail. For numerical studies we will use the central values of the following average PDG(20) numbers [35] for the running masses of c and b-quarks, namely m c (m 2 c ) = 1.27 ± 0.02 GeV, m b (m 2 b ) = 4.18 +0.03 −0.02 GeV. In accordance with the results of [91] obtained from the LHC tt experimental data and given in [39], for t-quark we assume m t (m 2 t ) = 164.3 ± 0.6 GeV that does not contradict the values of the running t-quark mass presented in PDG (20).
As the initial normalization point we take are defined using the corresponding N 3 LO matching transformation conditions, derived in [92,93] 16 (the corresponding NNLO conditions were obtained in [94,95] and are naturally taken into account by us), where the matching scales are fixed by the values of the MS-scheme masses presented above. Using the corresponding inverse logarithmic N 3 LO approximation for α s , we find: 16 The inclusion in the numerical analysis of the five-loop threshold effects investigated in [92,96,97] and of the five-loop contribution to the QCD β-function [61,62] does not affect essentially the numerical values of the pole masses of heavy quarks.
Note, that these numerical expressions are in agreement with those demonstrated in [35].
Taking into account the values given above, the exact results (11a), (12a-12c) and data from The terms in braces are the numerical contributions of the fifth and sixth orders gotten within the various estimate approaches considered by us in this work. Despite the fact that these values are approximate, they reflect the specific behavior of the high-order PT corrections to the relation between the pole and MS-scheme running masses of heavy quarks, viz • For the case of the c-quark the numerical PT QCD corrections in (39a) (which are starting to increase from the O(a 3 s ) level) are keeping on their asymptotic growth at the O(a 5 s ) and O(a 6 s ) orders. Indeed, the five-loop contribution is almost 2 times larger than the four-loop expression and the six-loop correction is more than 2 times greater than the five-loop one being estimated and is even larger than the first term of this PT series. This effect is strongly related to both the influence of the moderately large value of the coupling constant α s (m 2 c ) and to the renormalon contributions to the ratio M c /m c (m 2 c ). In this regard, in the modern high-precisely phenomenologically-oriented studies it is more appropriate to use the concept of the running c-quark mass that does not suffer from the factorial renormalon behavior.
• The estimates made for the b-quark signals that the asymptotic nature of Eq. (39b) is starting to reveal itself from the O(a 4 s )-contribution (except for the FL, M → m results where the O(a 5 s ) contribution is less than the O(a 4 s ) one). Note also that in the case of the application of the ECH-approach we observe the peculiar stabilization feature of the four-, five-and six-loop corrections to the on-shell-MS mass: within the accuracy considered by us in (39b) they coincide. However, due to the existing theoretical uncertainties discussed in Sect. 4 we do not consider this observation as a physical one.
• The relation (39c) demonstrates a decrease of the O(a 5 s ) and O(a 6 s )-corrections in all estimated approaches being investigated by us. Therefore, the asymptotic behavior of the PT series for the ratio M t /m t (m 2 t ) is not yet manifesting itself at the six-loop level.
One should emphasize that in fact the coefficients t M k depends substantially on a choice of the scale parameter μ (see Eqs. (2)(3) and [14,16,41,49] as well). For instance, shifting it from μ = m c to μ = 3 GeV one can delay the order of manifestation of the renormalon factorial growth in corrections to the ratio M c /m c (μ 2 ) (see e.g. [17,[20][21][22]24]) and move it to the fourth order of PT (as in the case of the b-quark).
Note also that the effects of the massive lighter quarks in the coefficients t M k are no less essential than the RGcontrollable ones responsible for the shift of the renormalization scale. They are rather important in both theoretical and phenomenological studies related to the determination of the charm, bottom and top-quark masses (see e.g. [22,25,27,28,89]). These massive effects were exactly calculated in [2,98] (see also the recent work [99]). However, in this paper we do not study the extra theoretical uncertainties related to the incorporation in analysis of the effects of massive lighter quarks.
Using the known results of [8] one may analyze the asymptotic structure of the relation between the pole and running t-quark masses in more details. Combining the six-loop FLexpression from (39c) with the results of [8] normalized at μ 2 = m 2 q and utilizing the NNA approximation, one can arrive to the following numerical representation for the topquark pole mass: This approximate expression indicates that in the case of the top quark the first traces of the asymptotic nature of the corresponding perturbative QCD series is observed above the seventh order of PT. Indeed, the contributions of the seventh, eighth and ninth orders are comparable to each other. Further, using the higher-order t-quark estimates for the coefficients t M k , obtained in [50] with help of the IRR-based formula (33), we conclude that in this case the numerical contributions to the relation between M t and m t (m 2 t ) are close to the FLones presented in (40). Therefore, the statement about the manifestation of the asymptotic behavior of the PT series for the ratio M t /m t (m 2 t ) after the seventh order seems to us quite reliable.

Conclusion
In this work we have estimated the values of the O(a 5 s ) and O(a 6 s )-contributions to the relation between the pole and running masses of heavy quarks. For this aim we have utilized three different approaches, namely the effective charges motivated method in its two variants (with and without the π 2 -effects of the analytic continuation from the Euclidean to Minkowskian region), the Naive-Nonabelianization procedure applied to the results of calculations of the leading renormalon-type terms and the technique based on the application of the renormalon asymptotic formula with the normalization factor fixed in two ways. The "kinematic" analytical continuation effects, which were modeled with help of the Källen-Lehmann type dispersion relation, have been associated by us with the π 2n -terms typical to the Minkowskian on-shell scheme.
As a result of these estimates we have obtained that at the fixed number of massless quarks the approximate values of the coefficients t M 5 and t M 6 evaluated by all three approaches are qualitatively consistent with each other (on average with a factor two). Further, using these results we have studied the flavor dependencies of these terms and established their signalternating character in n l (as in the case of the already exactly calculated two, three and four-loop ones). Herewith, we especially emphasize that upon studying of the n l -structure of the coefficients t M 5 and t M 6 estimated with help of the renormalon asymptotic formula the detailed information on the normalization factor N m in the expression (33) plays the essential role.
Further we have considered the asymptotic structure of the relation between the pole M q and MS-scheme running masses m q (m 2 q ) of the real heavy quarks. In comparison with the cases of the charm and bottom quarks where the asymptotic behavior manifests itself in the third and fourth orders of PT correspondingly, the asymptotic nature of the analogous PT series for the top quark seems to reveal itself above the seventh order. Therefore, the concept of the top-quark pole mass may be used safely in the modern phenomenologically and theoretically oriented studies.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Application of the least squares method
Let us discuss in more details the features of application of the least squares method. Following the studies done in [18] we use the results of semi-analytical calculations of the term t M 4 at the fixed number of massless quarks [17] and obtain the following overdetermined system of linear equations with two unknown parameters t M 4,0 and t M 4,1 normalized at the point μ 2 = m 2 q and defined in (5): Herewith, we restrict ourselves by the consideration of n l from the range 3 ≤ n l ≤ 15, where the lower bound is fixed by us keeping in mind that we analyze the behavior of perturbative series for the relation between the pole and running masses of heavy quarks, while the upper bound is following from the Banks-Zaks ansatz n l < 31/2 [100], which insures that in the considered region of n l the QCD asymptotic freedom property is not violated.
To apply the LSM for solving the system (A.1) one should first introduce the Φ-function, which is equal to the sum of the squares of the deviations of all equations in this system: (A. 6) In these expressions the second term under the square root reflects the effect of the correlation between uncertainties shown in the system (A.1).
However, if there were even no errors in (A.1), the LSM would still provide uncertainties related to the quality of the reproducing of the input data. These uncertainties can be directly calculated from the following formulas: Although the central values of the results (A.11a) and (A.11b) were obtained by different ways, they coincide. This is an additional argument in favor of the applicability of the least squares method and the reliability of the results of (A.11a).
where the uncertainties of the four-loop terms are the meansquare errors following from (B.12).
The expressions (B.13a-B.13c) demonstrate the asymptotic character of the corresponding perturbative QCD series. Indeed, all relations contain significantly growing and strictly sign-constant coefficients.
Turn now to the study of the PT series for the relation between the pole M l and running masses m l (μ 2 ) of the charged leptons (e, μ, τ ) in QED. In the case of the electron and muon their pole masses are the directly measurable parameters. In spite of the fact that the heavy τ -lepton is decaying rather fast, one can also introduce as its main characteristic the pole mass as well. It may be extracted, for instance, from the corresponding experimental data for the threshold behavior of the τ + τ − total cross-section production in e + e − collisions (see e.g. [107]). However, like in the case of quarks, it is also possible to define the MS-scheme running masses of the charged leptons which may be also used in the analysis of the experimental data [108]. Let us consider the structure of the ratio M l /m l (m 2 l ) in more details. Using the U (1)-limit of the results of the diagram calculations [1,2,[4][5][6][7] and [17]  where a = α(m 2 l )/π is the QED coupling constant defined in the MS-scheme and N l is the number of the massless charged leptons. The first four-loop N l -independent (the N ldependent one as well) term in the curly braces corresponds to the abelian U (1)-limit of the results of the semi-analytical calculations [17] carried out for the case of the SU (N c ) theory with n l massless quarks. The second one follows from the recent high-precision (about 1100 digits) four-loop computations of the on-shell mass renormalization constant Z OS m performed in QED in [72] (see also [109] where the wave function renormalization constant Z OS 2 was also found) for the case of N l = 0. It is worth emphasizing that the results of [72] are in very good agreement with the ones following from [17].
Taking into account Eq. (B.14) and keeping in mind that for the cases of the electron, muon and τ -lepton one should set N l = 0, 1, 2 respectively, we arrive to the following expressions: where the uncertainties of the four-loop terms are the meansquare errors following from (B.14). Unlike Eqs. (B.13a-B.13c) the expressions (B.15a-B.15c) demonstrate the absence of any regular sign-constant or signalternating structure of the related PT series (besides the case of the τ -lepton). The same feature is observed when the running MS-scheme QED parameters (the masses of charged leptons and coupling constant) are normalized at the scale μ 2 = M 2 l (in this case the sign-alternating structure of the relation between M τ and m τ (M 2 τ ) is not manifested itself anymore). Therefore, the point of view appearing from time-to-time in the literature that the asymptotic perturbative series in the QED should have sign-alternating structure, which is based in part on the theoretical studies presented in [102,103], seems to be not a general rule. Note that the sign-alternating behavior is realized nowadays only in the perturbative expression for the anomalous magnetic moment of electron, which is known at present with high precision up to the five-loop term (for the most recent results of its numerical evaluation see [110,111]).
One should mention that the issue of the non-regular behavior of the corrections to the relation between the pole and running masses of the charged leptons in QED was first raised in [112] upon the three-loop analysis of the numerical expressions for this relation being analytically evaluated later on in [113]. It may be of interest to understand better the discussed non-regular structure of the asymptotic QED series in the future.