Renormalon-motivated evaluation of QCD observables

A method of evaluation of spacelike QCD observables ${\cal D}(Q^2)$ is developed, motivated by the renormalon structure of these quantities. A related auxiliary quantity ${\widetilde {\cal D}}(Q^2)$ is introduced, which is renomalization scale independent only at the one-loop level, and a large-$\beta_0$-type renormalon motivated ansatz is made for the Borel transform of this quantity. This leads to a `dressed' Borel transform of the considered observable ${\cal D}(Q^2)$. From there, a Neubert-type distribution is obtained for the observable. The described method is then applied to the massless Adler function and the related decay ratio of the $\tau$ lepton semihadronic decays. Comparisons are then made with an evaluation method at higher truncated orders, developed in our earlier works, which is a renormalization scale invariant extension of the diagonal Pad\'e approximants.


I. INTRODUCTION
The perturbative QCD (pQCD) running coupling a(Q 2 ) ≡ α s (Q 2 )/π (where Q 2 ≡ −q 2 ), in the usual theoretical known renormalization schemes, has the peculiar property of having a significantly different regime of holomorphic (analytic) behavior in the Q 2 -complex plane than the spacelike QCD observables D(Q 2 ) such as current correlators, nucleon structure functions and their sum rules. Namely, on the one hand, the general principles of Quantum Field Theories (QFTs) imply [1,2] that these observables D(Q 2 ) are holomorphic functions in the Q 2 -complex plane with the exception of a part of the negative semiaxis, Q 2 ∈ C\(−∞, −M 2 thr ], where M thr ∼ 0.1 GeV is a threshold scale comparable with the lightest meson mass. 1 On the other hand, the pQCD coupling a(Q 2 ) has in general singularities along the negative semiaxis and, in addition, singularities outside the negative semiaxis, usually on the positive semiaxis 0 ≤ Q 2 ≤ Λ 2 Lan. . These are called Landau singuarities or Landau ghosts, and the point Q 2 br. = Λ 2 Lan. ∼ 10 −2 -1 GeV 2 is the Landau branching point. This difference between the singularities of D(Q 2 ) and a(Q 2 ) represents a problem in QCD from the theoretical, and from the practical point of view. Theoretically, in pQCD the leading-twist part (and the Wilson coefficients of the higher-twist part) of D(Q 2 ) is a function of a(κQ 2 ), where 0 < κ ∼ 1 is a chosen fixed renormalization scale parameter. Hence the evaluated total observable D(a(κQ 2 ); Q 2 ) ev. does not fulfill the holomorphic properties required by the QFT principles. Practically, for low spacelike scales Q 2 , |Q 2 | 1 GeV 2 , the argument µ 2 = κQ 2 in the coupling a(κQ 2 ) is either close to or within the regime of the Landau singularities, making the evaluation of a(κQ 2 ) and thus of D(a(κQ 2 ); Q 2 ) ev. either unreliable or impossible.
Other holomorphic variants of QCD have been proposed and used since then, cf. Refs. . 2 The significant difference between most of them and the (F)APT is usually the behavior of the spectral function of the coupling ρ A (σ) in the low-momentum regime 0 ≤ σ 1 GeV 2 . The coupling of Refs. [15] is infinite at Q 2 = 0, but in most of the other works the coupling is finite nonzero there, 0 < A(0) < ∞. Nonetheless, in some of the works the coupling has zero value at Q 2 = 0, A(0) = 0, cf. Refs. [17,29,34,35,40]. On the other hand, large volume lattice calculations of the dressing functions of the Landau gauge gluon and ghost propagators [41][42][43][44][45] give at low Q 2 the results which can be interpreted as implying A(0) = 0, when the coupling is defined in a natural way as the product of the obtained gluon and the square of the ghost dressing function, A(Q 2 ) = const × Z gl (Q 2 )Z gh (Q 2 ) 2 . Similar behavior of the mentioned dressing functions is also obtained in the Gribov-related and Dyson-Schwinger Equations (DSE) approaches [45][46][47]. A different definition [48,49] of the low-Q 2 running coupling A(Q 2 ), which involves an additional factor Q 2 /[Q 2 + M (Q 2 ) 2 ], where the parameter function M (Q 2 ) is usually called dynamical gluon mass, leads with such dressing functions to a coupling with A(0) > 0.
In most of the mentioned QCD variants (AQCD), the coupling A(Q 2 ) and its power analogs A n (Q 2 ) are, or can be, constructed by dispersive methods, automatically ensuring the wanted holomorphic properties of the coupling. In such dispersive approaches, nonperturbative contributions are naturally generated in the couplings A(Q 2 ) and A n (Q 2 ). Similar kind of dispersive approaches can be applied also directly to the spacelike QCD observables D(Q 2 ) [33,40,[50][51][52][53][54][55][56], i.e., without referring to the coupling.
In this work, a method motivated by the renormalon structure of spacelike QCD observables D(Q 2 ) is developed. In order to obtain a practical and unambiguous evaluation of the observables D(Q 2 ) with such a method, the coupling A(Q 2 ) should not have Landau singularities. As the starting point, an auxiliary quantity D(Q 2 ) is introduced in Sec. II, which is renormalization scale invariant only at the one-loop level. Physically motivated one-loop-type (largeβ 0 -type) renormalon ansätze are made for the Borel transform B[ D](u) of the latter quantity, leading to "dressed" renormalon expression for the Borel transform B[D](u) of the original observable. In Sec. III, the parameters in the expression for B [ D](u) are fixed by requiring the reproduction of the known perturbation series coefficients of the observable D(Q 2 ), specifically in the case of the Adler function. Subsequent application of the Neubert approach [57] to the Borel transform of this observable gives us the characteristic distribution function G D (t) of this observable. With G D (t) we can evaluate the value of D(Q 2 ), the evaluation being without infrared (IR) ambiguity in AQCD variants with IR-safe couplings A(Q 2 ). In the usual perturbative QCD (pQCD) the evaluation has ambiguity due the concurrence of the Landau singularities of the pQCD coupling a(Q 2 ) = α s (Q 2 )/π and the IR renormalon. This approach is then applied to the evaluation of the massless Adler function D(Q 2 ) in Sec. III B and the related (timelike) τ lepton semihadronic decay ratio r τ in Sec. III C, in two versions of QCD with IR-safe coupling (2δ [23,24] and 3δ AQCD [35]), and in pQCD in the corresponding schemes. In Sec. IV the obtained results are compared with those obtained with a generalization of the diagonal Padé method at high truncated orders, a method developed in our earlier works. The conclusions are presented in Sec. V. Some additional details are presented in Appendices: in Appendix A specific recursion relations for the coefficients k m (n) and k m (n) introduced in Sec. II A; in Appendix B the two used AQCD variants; in Appendix C the method of obtaining specific constants C (X) i,j appearing in Sec. II C is explained; in Appendix D explicit expressions are given for specific integrals of the Adler characteristic functions needed in Sec. III C.
Repeated use of the RGE (3) makes it possible to express the logarithmic derivatives Eq. (2) as a power series a n (Q 2 ) = a(Q 2 ) n + ∞ m=1 k m (n) a(Q 2 ) n+m .
Step-by-step inversion of these relations makes it possible to express separate powers in terms of the logarithmic derivatives a(Q 2 ) n = a n (Q 2 ) + ∞ m=1 k m (n) a n+m (Q 2 ) .
For example, we have etc. When these expressions for each of the powers a n are substituted in the power expansion Eq. (1) of D(Q 2 ), the rearranged expansion in the logarithmic detivatives (logarithmic perturbation expansion -'lpt') is obtained where the new coefficients d n are specific combinations of d n , d n−1 ,. . ., d 0 and where the coefficients k s (n + 1 − s) are those appearing in the relations (6). The more explicit form of Eqs. (9) is where the coefficients k s (n + 1 − s) are those appearing in the relations (4). We refer to Appendix A for recursion relations which allow us to obtain the coefficients k m (n) and k m (n) to any order in any given renormalization scheme.
Having obtained the modified coefficients d n [via Eqs. (9)-(10)], an auxiliary quantity D(Q 2 ) is constructed whose power expansion is obtained by the formal replacements a n (Q 2 ) → a(Q 2 ) n in the logarithmic perturbation expansion (8) While the observable D(Q 2 ) is independent of the renormalization scale µ 2 (i.e., κ-independent where κ ≡ µ 2 /Q 2 ) the quantity D is not an observable because it is κ-independent only at the one-loop level. 3,4 Namely, the dependence of the coefficients d n (κ) and thus of d n (κ) on the renormalization scale parameter κ is determined uniquely by the µ 2 -independence of D pt (Q 2 ), and then the quantity is κ-dependent at the level beyond one-loop. More specifically, it is possible to check that µ 2 -independence of the series (13b) implies the following κ-dependence of the coefficients d n (κ): As a consequence, d 0 (κ) is κ-independent, and where d n ≡ d n (1). This implies that, if a(µ 2 ) were obtained from a(Q 2 ) by one-loop RGE running, the expansion (14) would be µ 2 -independent (note: κ ≡ µ 2 /Q 2 ). Based on the relations (15), it is the straightforward to show that the Borel transform of the quantity D(Q 2 ; κ) of Eq. (14) has a simple ("one-loop-type") κ-dependence We point out that this is not an approximate, but exact κ-dependence. On the other hand, the Borel transform of the perturbation series Eq. (13a) of the observable D(Q 2 ) has the κ-dependence that is only approximately of the type (18), namely at the one-loop level approximation. 3 In this work the κ-dependent quantities at κ = 1 (i.e., for µ 2 = Q 2 ) are simply denoted without reference to κ; e.g., dn(κ = 1) ≡ dn, . 4 It is straightforward to see that the relations (9) and (11) are valid at any κ (> 0), not just κ = 1, i.e., the coefficients ks(m) and ks(m) are independent of κ: dn(κ) = n−1 s=0 ks(n + 1 − s) d n−s (κ) and dn(κ) = n−1 s=0 ks(n + 1 − s) d n−s (κ).
These considerations suggest that the Borel transform B[ D](u) of the auxiliary quantity D(Q 2 ) has a one-loop (large-β 0 ) type renormalon structures where p and k are positive integers. The first structure represents a large-β 0 -type k-fold IR renormalon at u = p, and the second a large-β 0 -type k-fold UV renormalon at u = −p. For example, the IR and UV p double pole renormalons (DP, k = 2) generate coefficients d n ∝ (n + 1)!(β 0 /p) n and (n + 1)!(−β 0 /p) n , respectively; the single pole (SP, k = 1) variant generates d n ∝ (n)!(β 0 /p) n and (n)!(−β 0 /p) n , respectively. In addition to these structures, the IR and UV p subleading renormalon (SL, 'k = 0') can (and will) be included which, although not meromorphic functions, generate the analogous "subleading" type of coefficients d n = (n − 1)!(β 0 /p) n and (n − 1)!(−β 0 /p) n , respectively. In practice, for QCD spacelike observables D(Q 2 ) we know the exact values of only the first three or four coefficients d n and d n (i.e., n = 0, 1, 2, 3). One of the aims of the present work is to obtain a physically motivated estimate of all the other higher-order coefficients d n and d n (n ≥ 4). The idea is to make a physically motivated ansatz for the Borel transform B[ D](u) as a sum of the (one-loop type) renormalon terms of the form Eqs. (19)- (20), and adjust the weights of these renormalon terms so that the (known) values of the first three or four coefficients d n are reproduced. 5 However, before embarking on this, it has to be checked first whether the forms Eqs. (19)- (20) generate physically acceptable coefficients d n . Namely, the generated coefficients d n lead, via the relations (11), to the coefficients d n for the power expansion (1) of the observable D(Q 2 ), 6 and the question is whether the (higher-order) coefficients d n obtained in this way fulfill the physically motivated expectations. The latter expectations are contained in the expected renormalon structure of the Borel transform of the observable D(Q 2 ) B. Full renormalon structure of D(Q 2 ) In many cases of spacelike observables D(Q 2 ) such as Adler function, the expected structure of B[D](u; κ) is known and comes from the IR and UV renormalons. The IR renormalon leads to an ambiguity of the Borel integral, the origin of this ambiguity is the low-momentum regime of the Feynman integrals, and consequently [58][59][60] the Q 2dependence of these terms must agree with the power-suppressed (higher-twist) terms appearing in the Operator Product Expansion (OPE). A D-dimensional OPE term of the spacelike observable D(Q 2 ) has the form where which then gives the relation between O D (Q 2 ) and Ô D Furthermore, RGE gives for the inverse powers of Q 2 the expression where the overall constant is Q 2 -independent and renormalization scheme independent, and the b The Q 2 -dependence of the OPE operator term is then where the first series, with coefficientsĉ . .) and of the subleading effects of the exponent of the integral with the anomalous dimension (24) On the other hand, the Q 2 -dependence coming from the IR renormalon ambiguity of the IR u = p renormalon can be obtained by evaluating the imaginary part ImD(Q 2 ) (±) IR,p,BI of the generalized Principal Value of the Borel integral where the Borel transform of the IR p renormalon takes the ansatz Direct evaluation gives for this ambiguity (cf., e.g., [61]) As mentioned, the Q 2 -dependence of the OPE term Eq. (27) is the same as that of the renormalon ambiguity Eq. (31), leading to the following expressions of the Borel renormalon parameters in terms of the corresponding OPE operator parameters: ).
The coefficient p is a positive integer for spacelike observables. The corresponding expressions for the coefficients (d n ) IR,p at a n+1 Q in the perturbation expansion ('pt') can then be obtained directly from the Borel transform (30), using the expansion (21) for κ = 1 When considering the UV u = −p renormalon contribution D(Q 2 ) UV,p , whose Borel transform has the ansatz the previous analysis can be repeated analogously [62] (cf. also [63]), considering a Q as negative, a Q = −|a Q |, and in the Borel integral the integration over u now goes from ±i to ±i − ∞. The corresponding (formal) OPE operators have dimension D = −2p (< 0), and the relations analogous to those of Eqs. (32) are and the corresponding coefficients (d n ) UV,p at a n+1 Q in the perturbation expansion ('pt') of D UV,p are C. Generation of the full renormalon structure from the one-loop-type renormalon structure This Section addresses the question of what kind of the renormalon structure for the observable D(Q 2 ) is obtained when the one-loop-type of renormalon structures Eqs. (19)- (20) are adopted for the auxiliary quantity D(Q 2 ) Eq. (12). We will see that the full renormalon structures of the type Eqs. (30) with (32) are obtained in the case of IR renormalons, and those of Eqs. (35) with (36) in the case of UV renormalons. The numerical investigation will be performed in more detail in two specific renormalization schemes: the Lambert scheme with c 2 = −4.9, and (four-loop) Lambert MiniMOM (LMM) scheme. These are the two schemes in which IR-safe (and holomorphic) QCD couplings A(Q 2 ) were constructed [23,35] which can be used naturally later on in the numerical evaluations of the observables D(Q 2 ), giving unambiguous numerical results. This is so because these couplings, although coinciding practically with the underlying perturbative coupling a(Q 2 ) at high |Q 2 |, do not have Landau singularities at lower |Q 2 | 1 GeV 2 . For more details on these couplings, the reader is referred to Appendix B.
In practice, the Borel transforms of the form Eqs. (19)-(20) generate d n coefficients from which, in practice, the coefficients d n up to n = n max = 70 were generated via the relations (11), using Mathematica software [64] (cf. also Appendix A). These d n coefficients, at large n, are then compared with the expressions of the form of Eqs. (34) and (37) originating from the form of the full renormalon Borel transforms B[D](u) as suggested by the theory, Eqs. (30)- (32) and (35)- (36). It will be seen that the numerical results confirm that these structures are really reproduced.
These Borel transforms generate the following perturbation expansion coefficients d n of the corresponding auxiliary quantities D(Q 2 ): The corresponding coefficients d n of the observable D(Q 2 ) that are obtained numerically from these d n 's by the aforementioned relations (11), turn out to agree at high n with the following expressions to a high precision (n ≤ n max with n max = 70 was used): where the coefficients b (D) j are given in Eqs. (26), and for the index γ p the following notation is used [cf. also Eq. (32a)]: The coefficients C   2,k , and the ratios d X p,k / d X p,k , for X=IR with p = 2, 3 (D = 2p), and X=UV and p = 1 (D = −2p) in two renormalization schemes: (four-loop) Lambert MiniMOM (LMM) scheme, and Lambert scheme (Lamb., with c2 = −4.9).
type LMM : C The Borel transforms corresponding to the expressions (40) are O2p /β 0 in Eq. (32a) is an integer. The latter condition appears to be satisfied for spacelike observables.
Analogously, also the UV p = 1 renormalons were considered numerically, in the mentioned three renormalization schemes, for the double-pole (DP), single-pole (SP), and the subleading (SL) cases: They generate the coefficients respectively, 7 for d UV 1,j = 1/π. From these coefficients d n , the coefficients d n generated by the relations (11) agree  Table I, but now for the LMM scheme truncated at c4 (TLMM: cj = 0 for j ≥ 5), and for the MS scheme (cj = 0 for j ≥ 5). In both schemes N f = 3 was taken.
type TLMM : C numerically with the following expressions to a high precision: and for the index γ p (with p = 1) the following notation was used [cf. Eq. (36a)]: The generated Borel transforms are The numerically determined coefficients C (−2) 1,j and ratios d UV 1,j / d UV 1,j are given in the lower part of Table I (UV, p = 1), in the two mentioned renormalization schemes.
One can also ask how the values of the coefficients and ratios in Table I would be affected if the renormalization scheme were truncated, e.g., c j = 0 for j ≥ 5. The results for the LMM scheme thus truncated are presented in Table  II, where also the results are included for the MS scheme (with N f = 3) which is known to c 4 . Comparing the results for LMM in Table I and for its truncated version in Table II, one can see that the truncation does not appreciably affect them. The UV p = 1 results are almost unaffected by this truncation.
The results in Tables I and II indicate that, in the cases of SP (simple pole of B[ D](u)), the correction coefficients C (D) j,1 are close to zero or even compatible with zero. Furthermore, the ratios d X p,k / d X p,k in each of the considered cases  Table III shows, in three different renormalization schemes (LMM, c 2 = −4.9 Lambert, and MS), the convergence of the generated d n coefficients in the case of IR p = 2 SP, with increasing n. Specifically, d n coefficients are those given in Eq. (39b) with p = 2, the coefficients d n are then generated via the relations (11) up to 8 n = 70, and they are divided by the truncated versions of the n-dependent part J(n) of Eq. (40b) with p = 2 The values of the coefficients C  Tables I and II. The results presented in Table III show strong convergence of the mentioned ratios, especially d n /J(n) (NNLO) , toward n-independent values when n increases.

A. Construction of the generating Borel transforms
The massless Adler function, which is a logarithmic derivative of the light-quark current correlator, is a specific example of a spacelike QCD observable for which we have a large amount of theoretical information available. Namely, its perturbation expansion (1) is known up to ∼ a 4 , i.e., the coefficients d n for n ≤ 3 are exactly known [65][66][67]. Further, the large-β 0 (LB) expansion of its Borel transform is also known [59,68,69], showing 9 that in the large-β 0 (resummed one-loop) approximation the IR renormalon poles are double for p = 3, 4, . . .
]. In addition, for the IR p = 2 renormalon pole the subleading coefficientĉ Using all this information, a physically motivated ansatz for the Borel transform B[ D](u) of the auxiliary quantity D(Q 2 ) Eq. (12) of the Adler function will be written, where the IR p = 2 leading and subleading renormalons are included, as well as the IR p = 3 and UV p = 1 leading renormalons Here, the superscript '(4 P)' indicates that the ansatz contains four adjustable parameters: the "scaling" parameter K and the renormalon residue parameters d IR 2,1 , d IR 3,2 and d UV 1,2 . It turns out that the IR p = 2 subleading parameter α is fixed by the knowledge of the subleading coefficientĉ Eq. (50). Each one-loop-type IR renormalon term in the ansatz (51) generates the Borel transforms of the type (42) and the corresponding contributions (40) to the perturbation coefficients d n . When assuming K = 0 and requiring that the two IR p = 2 renormalons in Eq. (51) together generate the (d n ) IR,p=2 coefficient of the form Eq. (34) with the subleading part there havingĉ (4) 1 as given in Eq. (50), the following condition is obtained: where .
The ratios and the parameter C 1,1 appearing in these expressions are given in Table I. On the basis of the results of the previous Section, the expression (51) generates the following Borel transform of the full Adler function D, at the 9 The MS scale convention is used throughout this work. We recall that dn(κ) can be expanded in powers of β 0 , namely dn(κ) = cn,n(κ)β n 0 + c n,n−1 (κ)β n−1 0 + . . . + c n,−1 (κ)β −1 0 , and d (LB) n (κ) = cn,n(κ)β n 0 , where cn,n(κ) is renormalization scale dependent and renormalization scheme independent [21,70], i.e., independent of the scheme parameters c j = β j /β 0 (j ≥ 2). renormalization scale µ 2 = κQ 2 = exp(− K) Q 2 : and the notations Eqs. (41) and (46) for γ p and γ p were used. The above Borel transform generates the Adler function coefficients d n (κ) at the value of the renormalization scale parameter κ (≡ µ 2 /Q 2 ) = exp(− K). As mentioned, if the value of the renormalization scale parameter were κ = 1 ( K = 0), then the requirement of the reproduction of the correct contribution (d n ) IR,p=2 of Eq. (34), at subleading order, with the known value ofĉ 1 Eq. (50), would imply that the α coefficient at d IR 2,1 in Eq. (54) must have the value as given in Eq. (53). It turns out that the effect of κ = 1 ( K = 0) does not change this relation and the relation (52). This is so because i.e., the change of renormalization scale changes the coefficient only by a constant factor, with no subleading corrections (but with subsubleading corrections). This relation can be understood if the IR p = 2 part of the Borel B[ D](u) of Eq. (51) is reexpressed by expanding the exponential exp( Ku) around u = 2, which gives The term (1 − u/2) ln(1 − u/2) in the brackets is subsubleading (SSL); comparison of the results Eqs. (40d) with (40b) then gives that this term gives relative corrections ∼ 1/n 2 , i.e., the relation (55). Hence, in retrospect, we see that the relations (52)- (53), which are subleading in their nature, are not affected by K = 0 which is a subsubleading effect [apart from the overall factor exp(2 K)].
In the case of the mentioned renormalization schemes LMM and Lambert c 2 = −4.9 schemes, applicable in the 3δ AQCD and 2δ AQCD, respectively, and the truncated (at c 4 ) TLMM and MS schemes, the relations (52)-(53) and the numerical results of Tables I and II give us the values By far the most dominant source of uncertainty of α is the uncertainty δC  Tables I and II). Since the Borel transform B[ D](u) Eq. (51) has four free parameters, and its power expansion generates the coefficients d j , the four free parameters can be determined by the knowledge of the first four coefficients d j (j = 0, 1, 2, 3) in the considered scheme, whose values are given in Table IV ( Table V (first line). This model for the Adler function will  be applied for the 3δ AQCD [35] because this QCD variant was constructed in the LMM renormalization scheme. An additional goal here is to make comparisons of the results obtained in this way with the analogous results obtained in 2δ AQCD [23], which is a QCD variant constructed in the Lambert c 2 = −4.9 renormalization scheme. Therefore, a model for the Adler function will now be constructed in this Lambert renormalization scheme, by requiring additionally that the same value d 4 (MS, N f = 3) = 338.2 be generated in this model. This now implies one more condition (in total five conditions) in this scheme. Hence, for the c 2 = −4.9 Lambert renormalization scheme, the following ansatz with one more parameter is written (the IR p = 3 renormalon will have the double and the single pole): Applying the five conditions, the five parameters of this ansatz ( K and the four renormalon residues) are determined uniquely, and this solution is given in the second line of Table V. The same procedure was repeated for the MS scheme (no comparable AQCD version is available in that scheme, though), and the results are included in Table V. 10 Now we have the two models of Adler function, in 4LMM and Lambert schemes, the two models being presumably comparable because they give the same value of the ∼ a 5 perturbation coefficient d 4 , namely d 4 (MS, N f = 3) = 338.2.

B. Construction of the characteristic distribution function GD(t)
Now the characteristic function F D (t) will be constructed, for the spacelike observables D with a rather generic form of the Borel transform B[ D](u) of the correponding auxiliary quantity D. In this approach, mainly Neubert's construction [57] will be followed, who constructed such characteristic functions in the framework of the large-β 0 approximation in pQCD.
The characteristic function F D (t) is a function of a dimensionless parameter t (t > 0) such that the integral represents the leading-twist part of the observable D(Q 2 ), in the sense that it generates the correct logarithmic perturbation expansion (13b) of the observable when the coupling a(tQ 2 ) is Taylor-expanded around a(Q 2 ) 10 The polynomial five-loop β function [76] was used in the N f = 3 regime, using as the N f = 3 reference valueā 0 ≡ a(Q 2 0 ; MS) N f =3 = 0.0846346 at Q 2 0 = (2mc) 2 (where mc = 1.27 GeV). See Appendix B for more details, in particular footnote 23 there for comparison of numerical values.
where the notation Eq. (2) is used. This means that F D (t) must satisfy the following string of relations: Using these relations (with κ = 1) and the expansion (17) for the Borel transform B[ D](u), one obtains where the integration is in the complex u-plane paralell to the imaginary axis, and u 0 is any real value where the integral (62) exists, i.e., in the case of the Adler function one can take −1 < u 0 < +2. One can choose u 0 = +1, and write the above integral along the real axis in terms of the variable z such that u = 1 − iz For the Borel transforms B[ D] of the one-loop-type form (38) and (43), the integrals (64) for F D (t) can be evaluated in a straightforward way, by performing integration along judicially chosen contours and using Cauchy theorem. For the rather generic case of the Borel transform the following result for the integrated observable D(Q 2 ) is obtained: where the two characteristic functions are where Θ is the Heaviside step function, i.e., the result (66) can be written as The characteristic function G D (t) (= F D (e − K t)) is obtained by closing the integration path of z with the large semicircle in the upper (lower) half plane for t < 1 (t > 1). For the subleading (SL) contribution, the function G (SL) D (t) is obtained by closing the z-path in the upper half plane and integrating there along both sides of the cut (+i, +i∞) in order not to enclose it. It can be further noted that the SL contribution contains in the integrand the subtraction [a(te − K Q 2 ) − a(e − K Q 2 )] instead of simply a(te − K Q 2 ), because the SL contribution starts at ∼ a 2 because there d 0 = 0 (the corresponding Borel transform expansion starts at ∼ u 1 ). The relations (61) can now be rewritten in terms of G D (t) (= F D (e − K t)) and include the SL contribution It should be pointed out that the coupling a(te − K Q 2 ) in the integrand (66) is, in principle, running to any chosen loop order, thus representing the (leading-twist of the) full observable D(Q 2 ). If the coupling a(te − K Q 2 ) is one-loop running then the integral (66)[⇔ (68)] reproduces the perturbation expansion (12) of the auxiliary quantity D(Q 2 ) when the coupling (70) in the integral is expanded in powers of a(Q 2 ). The integral (66) is in pQCD in general ambiguous for Q 2 > 0 because of the Landau singularities of the pQCD coupling a(te − K Q 2 ) at low t values. To avoid this ambiguity, the integral (for Q 2 > 0) is evaluated with the path scale slightly above the real positive axis, a(te − K Q 2 + i ) and taking the real part of it (i.e., the generalized Principal Value) The ambiguity is proportional to the imaginary part On the other hand, the situation is essentially different in the QCD variants with IR-safe coupling A(Q 2 ) [the analog of a(Q 2 )]. In those variants, the coupling A(Q 2 ) is in general holomorphic (analytic) function in the Q 2 -complex plane with the exclusion of (part of) the negative semiaxis: Q 2 ∈ C\(−∞, −M 2 thr ], where M thr ∼ 0.1 GeV is a threshold scale generally of the order of the light meson mass. Two representative cases are the 2δAQCD and 3δAQCD, cf. Appendix B for more details. In such QCD variants, the evaluation of the integrals (66) [⇔ (68)] is unambiguous. In such AQCD frameworks, these integrals represent an unambigous resummation of the leading-twist part of the spacelike observable D(Q 2 ), if one takes the position that the Borel transforms B[ D](u) of the type (65) represent the correct generators of all the d n perturbation (lpt) coefficients of D(Q 2 ) Eq. (13b). As argued in Appendix B, most of the relations in Sec. II A survive in these QCD variants, with the substitutions: a → A, a n → A n , and a n → A n . Since the coupling A(Q 2 ) in these AQCD frameworks differs from the (underlying) pQCD coupling a(Q 2 ) in the same renormalization scheme by nonperturbative contributions, we have one (crucial) difference, namely that the power analogs A n are not simple powers of A (A n = A n ), cf. Appendix B. For example, the series (13) for the Adler function (and any other spacelike observable) in such AQCD frameworks gets the form where d 0 = d 0 = 1, the definition of the couplings A k and A k is given in Appendix B [Eqs. (B11) and (B14)] and, as always, κ ≡ µ 2 /Q 2 is the (arbitrary) dimensionless renormalization scale parameter (0 < κ ∼ 1). The terms in the series (73) do not suffer from Landau singularities at low positive Q 2 (0 ≤ Q 2 1 GeV 2 ), in contrast to the terms in the pQCD series (13). However, the series (73) are asymptotically divergent, as are also the pQCD series (13). The corresponding resummation of the series D AQCD (Q 2 ), with the characteristic functions G D (t) and G In contrast to pQCD, these integrals are unambiguous because the coupling A(Q 2 ) has no Landau singularities, and there is no need to employ the (generalized) Principal Value approach Eqs. (71)- (72). In Figs. 1(a),(b) the resulting Adler function D(Q 2 ) Ares is presented for positive Q 2 , as a function of Q ≡ Q 2 , in 3δ AQCD (in the LMM renormalization scheme) and in 2δ AQCD (in the Lambert c 2 = −4.9 renormalization scheme), respectively. In both cases, the results in the corresponding underlying pQCD (i.e., pQCD in the LMM and the Lambert renormalization scheme) are included, using the resummation form Eq. (71) with the uncertainty estimate Eq. (72). The pQCD curves are not included for very low Q < 0.6 GeV because there they have more erratical behavior.  [23] and in 3δ AQCD [35] in specific approaches. 11 The two programs used for the numerical evaluation are freely available at the www site [38], and are written in Mathematica language. Further, as mentioned in Sec. III A, the higher order d n coefficients of the Adler function in the three cases (LMM, Lambert and MS) do not mutually agree when transformed to a common renormalization scheme (e.g., to MS); however, the construction was such that the first five such coefficients mutually agree ( d MS 0 , . . . , d MS 5 ), the first four being the known coefficients ( d MS j ⇒ d MS j , where j = 0, . . . , 3; d 0 = 1), and the fifth coefficient being mutually the same (d MS 4 = 338.19). One aspect in Fig. 2 that catches the attention is the significantly different behavior of the radiative Adler function at Q < 0.5 GeV in the two AQCD variants; this is a consequence of the fact that the 2δ AQCD coupling A(Q 2 ) freezes in the deep IR regime at the nonzero value A(0) ≈ 0.66, while the 3δ AQCD coupling goes to zero as A(Q 2 ) ∼ Q 2 when Q 2 → 0, where the latter property is suggested by the results of the large volume lattice calculations of the dressing functions of the Landau gauge gluon and ghost propagators [41] (cf. also [42][43][44]). The reader can refer to Appendix B for more details on the couplings A and A n in 2δ and 3δ AQCD.

C. Evaluation of the semihadronic τ decay ratio rτ
In the previous Section it was described how to evaluate the leading-twist part of the spacelike observables D(Q 2 ) from the Borel transforms B[ D](u) of the auxiliary quantity D(Q 2 ). This evaluation is unambiguous in AQCD where no Landau singularities are present in the coupling A(Q 2 ) [and thus in A n (Q 2 )]. If D(Q 2 ) is the (massless) Adler function, this evaluation then allows one to evaluate, again unambiguously, the τ lepton semihadronic decay ratio 11 Cf. the next Sec. III C for more explanation on r where the Adler function D(Q 2 ) has expansions of the form Eqs. (13) with d 0 = d 0 = 1. 12 Using the expression (68) for D leads to The (IR-safe and holomorphic) coupling A can be written in terms of the dispersion integral along its cut where ρ A is the discontinuity (spectral) function of A along its cut, ρ A (σ) = ImA(−σ − i ), and the integration starts generically at σ = 0 [in 2δAQCD and 3δAQCD ρ A is zero up to σ = M 2 thr = M 2 1 , cf. Eq. (B3)]. Substituting the dispersion integral (77) into Eq. (76) and exchanging the order of integration, gives where where m 2 ≡ e − K m 2 τ . This then leads to the following expression for r (D=0) τ only in terms of the spectral function ρ A (σ) and the characteristic functions G D (t) and G (SL) D (t) of the Adler function: These double integrals can be reduced further to single integrals by integration by parts in t where we recall that m 2 = e − K m 2 τ ; the functions F j , F (SL) j and the constants C Using the expressions (67)

IV. FROM ASYMPTOTICALLY DIVERGENT SERIES TO A CONVERGENT SEQUENCE
Here, another method of evaluation of the leading-twist part of spacelike observables D(Q 2 ) will be used, a method which grows increasingly effective when the number of known coefficients d n in the logarithmic perturbation expansion (13b) increases. This method was proposed in Ref. [77] for the case when the number of known coefficients is even ( d 0 , . . . , d 2M −1 ), and was modified in Ref. [78] to be applicable when the number of known coefficients is odd. Later this method was applied in AQCD variants where the A(Q 2 ) coupling is IR-safe and holomorphic outside the negative semiaxis in the complex Q 2 -plane [25,70,79]. It is an extension of the diagonal Padé (dPA) approach, where the latter gives a result which is renormalization scale independent at the one-loop level [80].
where a µ ≡ a(µ 2 ) and the 2M − 1 coefficients N j , D k are determined by the condition that the expansion of the above expression in powers of a µ reproduce the (known) coefficients d 1 (κ), . . . , d 2M −1 (κ). The ratio (84) can always be decomposed into a sum of simple fractions where a (1−l.) (µ 2 ) (κ j Q 2 ) is the value of the coupling a RGE-evolved from the scale µ 2 (≡ κQ 2 ) to κ j Q 2 by one-loop RGE and the complex constants α j and κ j = κ exp( u j /β 0 ) in Eq. (85) are independent of the renormalization scale µ 2 ≡ κQ 2 and of the physical scale Q 2 , [77]. Further, we have M j=1 α j = 1.
The diagonal Padé (85) is independent of the initial renormalization scale µ 2 in the approximation of one-loop, but it is not exactly µ 2 -independent. To extend this resummation in such a way as to get the exact µ 2 -independence, the crucial point was to realize that α j and κ j are µ 2 -independent constants. The one-loop RGE-evolved quantities a (1−l.) (κ j Q 2 ) are replaced by the exactly (n-loop) RGE-evolved couplings a(κ j Q 2 ), leading to This is now exactly µ 2 -independent, 14 and can be shown [77] to agree with the perturbation expansion of the full observable D(Q 2 ) to the order ∼ a 2M This resummation was applied in pQCD [81] with reduced success because some of the complex (or real) parameters κ j usually turn out to be small (|κ j | < 1), so that the scale (κ j Q 2 ) is in the vicinity or on the Landau singularities, making evaluation a(κ j Q 2 ) unrealistic or impossible. Later, in Refs. [25,70,79], it was realized that this resummation can be applied in the AQCD frameworks, i.e., with IR-safe and holomorphic A(Q 2 ) coupling, with the same values of α j and κ j where D(Q 2 ) AQCD is the expansion in logarithmic derivatives A n (µ 2 ) in AQCD, Eq. (73a). In contrast to pQCD, in the AQCD variants the small |κ j | values are not a problem, because there are no Landau singularities of A(Q 2 ) [and thus of A n (Q 2 )] in the Q 2 -complex plane. This approach was applied to the known four-term truncated series D [4] Adl (Q 2 ) for the Adler function ( d j known up to j max = 2M − 1 = 3, i.e., M = 2) for various AQCD variants in Refs. [25,35,70,79]; and in Refs. [70,79] to the large-β 0 (LB) part of the Adler function, D (LB) Adl (Q 2 ), which is known to all orders [cf. Eqs. (49)]. However, the AQCD series (73) is in general asymptotically divergent. 15 Specifically, the sequence {D is asymptotically divergent at any Q 2 , because of the renormalon-type behavior of the coefficients d n ∼ n! at large n. Further, the extended dPA G Adl (Q 2 ), applied in various AQCD variants, the results of Refs. [70,79] strongly indicate that this is not so, and that the mentioned sequence of extended dPA's is a convergent sequence, for all complex Q 2 for which A(Q 2 ) is holomorphic, i.e., in the entire spacelike regime Q 2 ∈ C\(−∞, −M 2 thr ]. We recall that in Sec. III A the (leading-twist) Adler function D Adl (Q 2 ) was constructed in two models, i.e., the Adler function (73a) whose coefficients d n are generated: (a) by the one-loop type Borel B[ D](u) of Eq. (51), in the LMM renormalization scheme, where 3δAQCD [35] is applied; (b) by B[ D](u) of Eq. (58), in the c 2 = −4.9 Lambert renormalization scheme, where 2δAQCD [23] is applied. In Table V the parameters of these two Borel transforms were given. In both Adler function models, the first four known coefficients d j (⇔ d j ) are reproduced (j = 0, 1, 2, 3), and, in addition, both models give the same value of the next coefficient, which in the MS scheme is d 4 (MS, N f = 3) = 338.2. Further, we recall that in both models the full leading-twist Adler function can be evaluated (resummed) exactly, by integrations involving the corresponding characteristic functions G D (t) and G   In Table VI Table VI the corresponding values of the truncated AQCD series (91) were included in parentheses (with κ = 1), which clearly show, as expected, that such series are asymptotically divergent, with divergence setting in at terms  AQCD (Q 2 ; κ = 1) Eq. (B12), for N = 1, 2, . . . ('logderivs'); this series can be written also in the form of Eq. (B13) (with κ = 1), where the power analogs A n (Q 2 ) are not equal to the naive powers A(Q 2 ) n . For additional comparison, in Figs. 6 were included the results of the truncated AQCD series where the terms d n−1 A n (Q 2 ) are replaced by the naive (and thus incorrect) power terms d n−1 A(Q 2 ) n ('naive powers'). 16,17 It can be seen in Figs. 6 that the extended dPA sequence is consistently convergent, it converges to the exact value, and that the sequence of truncated AQCD is asymptotically divergent, i.e., it approximately stabilizes at 3 ≤ N ≤ 6, and for 6 < N is starts diverging. The sequence of the truncated series in naive powers, on the other 16 In [70] it was argued that the naive power terms A(Q 2 ) n in AQCD variants in general bring spurious uncontrollable nonperturbative contributions, in contrast to the logarithmic derivatives An(Q 2 ) and the related power analogs An(Q 2 ); cf. also Eqs. (B13)-(B14) in Appendix B. 17 On the other hand, in Ref. [82] a perturbative (pQCD) coupling a(Q 2 ) = A(Q 2 ) was constructed which has no Landau singularities and reproduces the correct values of rτ . In this case, the construction (B14) gives An ≈ a n (the equality becoming increasingly better when the truncation index N increases), and the extended Padé (88) or equivalently (90) can be applied, due to the absence of the Landau singularities, and gives a convergent sequence, presumably converging to the exact value. However, it is unclear how to construct a renormalon-motivated model (58) for the higher order 'exact' coefficients dn of the Adler function in this case, because the pQCD scheme of this holomorphic pQCD coupling has the scheme coefficients cn = βn/β 0 growing so fast with increasing n that the renormalon structure is severely affected by the transformation into this scheme (with the exception of the p = 1 UV renormalon).
hand, does not show any clear stabilization, it appears to be more divergent.
The τ decay ratio parameter r (D=0) τ can also be evaluated with this method, by evaluating the Adler function D(m 2 τ exp(iθ)) in the integrand of the contour integral Eq. (75) with the described extended dPA method. If the latter method gives convergent sequence in principle for all complex Q 2 [including those not far from the negative semiaxis where A(Q 2 ) has a cut], then it is expected that the obtained sequence r  83), although more slowly. This is really the case, as the obtained results presented in Table VII show.

V. CONCLUSIONS
In QCD we often face the problem of evaluation of the spacelike renormalization scale and scheme independent quantities D(Q 2 ), such as current correlators, nucleon structure functions and their sum rules, etc. The evaluation of the leading-twist part (which is usually dominant and strongly influences the values of the extracted parameters of the OPE higher-twist terms), has at least three aspects making the evaluation difficult and imprecise: (a) at low values |Q 2 | 1 GeV 2 , the evaluation of the pQCD coupling a(κQ 2 ) ≡ α s (κQ 2 )/π (where 0 < κ ∼ 1 is the renormalization scale parameter) cannot be performed reliably because the scale µ 2 = κQ 2 is either close to, or in the regime of, the Landau singularitiers of the coupling; (b) the coefficients d n (κ) of the perturbation series of D(Q 2 ) are not known, with the exception of the first few; (c) even if we knew these coefficients, or had a reasonable physically motivated estimate for them, the resulting perturbation series d n (κ)a(κQ 2 ) n+1 would be asymptotically divergent, even in the high-momentum regime |Q 2 | > 1 GeV 2 , due to the renormalon growth of the coefficients d n (κ) ∼ n!. This work addresses all three issues, and brings new insights in particular to the aspects (b) and (c).
The solution of the problem (a) has been known for some time, and it consists of two parts: (a1) replacing the pQCD coupling a(κQ 2 ) by a coupling A(κQ 2 ) which has no Landau singularities, i.e., it is holomorphic (analytic) in the complex Q 2 -plane with the exception of (a part of) the negative semiaxis, thus reflecting qualitatively the correct holomorphic properties of the spacelike observables D(Q 2 ) in the Q 2 -complex plane; (a2) reorganizing the perturbation series D(Q 2 ) = d n (κ)a(κQ 2 ) n+1 into a series D(Q 2 ) = d n (κ) a n+1 (κQ 2 ) with logarithmic derivatives a n+1 [Eqs. (2), (13)], and replacing these pQCD logarithmic derivatives by the corresponding logarithmic derivatives A n+1 of A [cf. Eqs. (73) and (B11)]. Such QCD variants can be called (holomorphic) AQCD frameworks.
In this work a solution to the problems (b) and (c) was proposed. In Sec. II, for any spacelike observable D(Q 2 ), an auxiliary quantity D(Q 2 ; κ) was introduced. Its perturbation (power) series D(Q 2 ; κ) = d n (κ)a(κQ 2 ) n+1 , Eq. (14), agrees with the (logarithmic derivatives) perturbation series of the observable D(Q 2 ) = d n (κ) a n+1 (κQ 2 ) only in the one-loop approximation (in which a n+1 = a n+1 ). The Borel transform B[ D](u; κ) of this quantity was shown to have the one-loop-type renormalization scale dependence, and consequently physically motivated specific ansätze for B[ D](u; κ) were proposed which capture the known (one-loop-type, or large-β 0 -type) renormalon structure of the observable D(Q 2 ). The parameters in this Borel transform were adjusted in such a way that the first few known coefficients d n of the observable D(Q 2 ) were reproduced. Reexpansion of the obtained B[ D](u) in powers of u then generated a (physically motivated) set of coefficients d n at all n. Reorganizing backward the obtained logarithmic derivative expansion of D(Q 2 ) = d n (κ) a n+1 (κQ 2 ) into power expansion D(Q 2 ) = d n (κ)a(κQ 2 ) n+1 then generated the perturbation (power) expansion coefficients d n at all n. The generation of d n coefficients from d n coefficients, via the replacement a n+1 → a n+1 in D, can be viewed as a dressing procedure which incorporates an important part of the radiative beyond-one-loop corrections. 18 This was further motivated in Sec. II C where the effect of this dressing was shown to give the full renormalon structure for the resulting Borel transforms B[D](u). Specifically for the (massless) Adler function D(Q 2 ), the perturbation expansion coefficients d n (and d n ) were generated in Sec. III A in three different renormalization schemes. This addresses the mentioned problem (b). 18 We notice that the dn coefficients, although generated by the one-loop-type Borel transform B[ D](u), contain also the radiative correction effects which are in general beyond-one-loop. Only the leading-β 0 (LB) part d However, according to the mentioned aspect (c), the resulting perturbation series for D(Q 2 ), Eqs. (13) in pQCD and Eq. (73) in holomorphic AQCD frameworks, are asymptotically divergent, primarily because of the fast growth of the coefficients d n , d n ∼ n!. Therefore, the Neubert-type of characteristic distribution function G D (t) of the observable D(Q 2 ) was constructed in Sec. III B, as the inverse Mellin transform of B[ D](u), leading to the evaluation (resummation) of (the leading-twist part of) the observable D(Q 2 ) as an integral over t of the product of G D (t)/t and the coupling a(tk Q 2 ) or A(tk Q 2 ), Eqs. (68) and (74) [where k = exp(− K) > 0]. In pQCD, due to the Landau singularities of the coupling a(tk Q 2 ), this evaluation becomes ambiguous at Q 2 > 0, cf. Eqs. (71)- (72), while in the holomorphic AQCD frameworks there is no ambiguity, cf. Eq. (74). In Sec. III B the (massless) Adler function as a function of Q 2 was explicitly evaluated in this way, in two AQCD frameworks, namely 3δ [35] and 2δ AQCD [23], as well as in the corresponding underlying pQCD frameworks and in MS pQCD. This addresses the aforementioned problem (c).
In addition, the described formalism was extended to the evaluation of the timelike observables R(s), when the latter can be expressed as integral transformations of the corresponding spacelike observables D(Q 2 ). Thus, in Sec. III C the semihadronic τ lepton decay ratio r It should be pointed out that evaluations with the described methods can be applied also to other spacelike and timelike QCD observables, in the QCD variants with IR-safe (and holomorphic) coupling A(Q 2 ), and will presumably lead to comparably favorable results as they do for the Adler function and τ decay ratio.
Once knowing the coefficients k i (j), the recursion relations for k m+1 (n) in terms of k p (r) (p ≤ m) are The pQCD running coupling a(Q 2 ) ≡ α s (Q 2 )/π is a solution of the perturbative RGE (pRGE) Eq. (3), where the first two β-coefficients, β 0 = (1/4)(11−2N f /3) and β 1 = (1/16)(102−38N f /3), are universal, i.e., scheme independent, in mass independent schemes. The other coefficients c j = β j /β 0 (j ≥ 2) characterize in pQCD the renormalization scheme [83]. Stated differently, the form of the function β(a; c 2 , c 3 , . . .) can be regarded as the definition of the renormalization scheme. The momentum scale parameter Λ QCD is not considered here as a scheme parameter, but rather as the definition of the momentum (re)scaling, and a scaling change can be equivalently described as a change of the renormalization scale. In this work, the MS scaling definition (Λ 2 QCD = Λ 2 ) is used throughout.
When integrating the pRGE in a given or chosen renormalization scheme, the resulting pQCD running coupling a(Q 2 ) usually acquires singularities on the positive axis in the Q 2 -complex plane, 0 ≤ Q 2 Λ 2 QCD (∼ 0.01-1 GeV 2 ), in addition to the otherwise expected singularities on the negative axis. On the other hand, the general principles of Quantum Field Theories imply that the spacelike QCD observables D(Q 2 ), such as current correlators and nucleon structure functions and their sum rules, are holomorphic (analytic) functions in the Q 2 -complex plane with the exception of a part of the negative semiaxis, Q 2 ∈ C\(−∞, −M 2 thr ], where M thr ∼ 0.1 GeV [1,2]. The pQCD running coupling a(Q 2 ) therefore usually does not reflect qualitatively these properties, because of the mentioned singularities (cut and branching points) on the positive axis, 0 ≤ Q 2 ≤ Λ 2 Lan. . This aspect of a(Q 2 ) is considered unfortunate, especially if the coupling a(Q 2 ) [or a(µ 2 ) with µ 2 = κQ 2 ∼ Q 2 ] is to be used to evaluate D(Q 2 ) at low values |Q 2 | 1 GeV 2 . These singularities are called Landau singularities or Landau ghosts, and the point Q 2 = Λ 2 Lan. is usually called the Landau branching point. The application of the Cauchy theory to the integrand a(Q 2 )/(Q 2 − Q 2 ) in the Q 2 -complex plane leads then to the following dispersion integral representation of the pQCD coupling a(Q 2 ): where ρ a (σ) = Ima(Q 2 = −σ − i ) is called the discontinuity or spectral function of a along its cut. The simplest elimination of the disturbing Landau singularities is to eliminate in the above dispersion integral the integration part along the positive σ (negative Q 2 ), this leading to the minimal analytic (or Analytic Perturbation Theory -APT) coupling [3][4][5] This coupling has the cut threshold σ min (≡ M 2 thr ) = 0, and the difference between this coupling and the underlying pQCD coupling a(Q 2 ) (i.e., the pQCD coupling in the same renormalization scheme) at large |Q 2 | > 1 GeV 2 is appreciable, namely A (APT) (Q 2 ) − a(Q 2 ) ∼ (Λ 2 QCD /Q 2 ) N with the index N = 1. We do know that at low squared momenta (|Q 2 |, σ 1 GeV 2 ) the pQCD coupling a(Q 2 ) and its spectral function ρ a (σ) do not describe correctly the physics. Therefore, it appears natural to replace in the dispersion integral (B2) the pQCD spectral function ρ a (σ) at σ 1 GeV 2 by another, unknown spectral function ρ A (σ) = ρ a (σ). An efficient way of parametrizing this spectral function ρ A (σ) in the low-momentum regime σ 1 GeV 2 is to represent it as a linear combination of Dirac delta functions. This suggests then the following form for the spectral function of the coupling A(Q 2 ): where we expect 0 < M 2 1 < . . . < M 2 n < M 2 0 , and M 2 0 ∼ 1GeV 2 can be called the pQCD-onset scale. The corresponding coupling is now The n Dirac delta functions in the spectral function thus give ∆A IR (Q 2 ) which is a linear combination of n simple fractions ∼ 1/(Q 2 +M 2 j ), and this can be represented as a near diagonal Padé approximant ∆A IR (Q 2 ) = [n/n−1](Q 2 ). Such Padé approximants are known to approximate usually the holomorphic functions in the Q 2 -complex plane increasingly well when the index n increases.
In Refs. [23,24] and [34,35], such couplings were constructed, with two (n = 2) and three (n = 3) Dirac delta functions, respectively, in specific renormalization schemes of the underlying pQCD coupling a. The (2n + 1) parameters (M 2 j , R j , j = 1, . . . , n; and M 2 0 ) were then fixed by several physically motivated conditions. Four of these conditions were obtained by requiring that the A(Q 2 ) coupling at high |Q 2 | > 1 GeV 2 practically coincides with the underlying pQCD In addition, at moderate momenta |Q 2 | ∼ m 2 τ (∼ 1 GeV 2 ) the requirement was imposed that the well measured physics of the semihadronic τ lepton decays be reproduced correctly, i.e., that the (massless and strangeless) τ decay ratio r These values were imposed in a specific chosen evaluation procedure that was considered relatively reliable. 20 Finally, the two additional conditions needed in 3δ AQCD were at very low momenta, namely that A (3δ) (Q 2 ) ∼ Q 2 when Q 2 → 0, and that A (3δ) (Q 2 ) acquires at positive Q 2 a local maximum at about Q 2 ≈ 0.135 GeV 2 , in the Lambert MiniMOM (LMM) renormalization scheme. These two conditions are suggested by the large volume lattice calculations [41] (cf. also [42][43][44]) of the dressing functions Z gl (Q 2 ) and Z gh (Q 2 ) of the Landau gauge gluon and ghost propagators in the MiniMOM (MM) renormalization scheme [84,85], where the lattice coupling was defined naturally as the product of these dressing functions: The Lambert MM (LMM) scheme, in which the 3δ AQCD coupling was constructed, is just the MM scheme of the lattice calculations [with the corresponding RGE scheme coefficients c j (MM), j ≥ 2] but with the usual (MStype) scaling (i.e., Λ QCD = Λ). In Ref. [34] the A (3δ) coupling was constructed in a scheme which agrees with the perturbative LMM scheme up to three-loop level (i.e., up to c 2 ); in Ref. [35] the construction was refined to agree with the known LMM scheme up to the four-loop level (i.e., up to c 3 ), i.e., the level known at present [84] for the MM scheme. On the other hand, the 2δ AQCD coupling A (2δ) (Q 2 ) was constructed in the underlying Lambert scheme with c 2 = −4.9 [24].
The β-functions β(a) defining the two schemes (Lambert and four-loop LMM), and thus the underlying pQCD 19 It is known that the higher dimensional (higher-twist) terms D > 0 contribute only little to rτ , cf. [35] and references therein. 20 In 3δ AQCD, which was constructed in the LMM renormalization scheme, r (D=0) τ was evaluated by the contour integration (75) with the Adler function evaluated with the truncated series in logarithmic derivatives, Eq. (73a), truncated at the last fully known coefficient [ d 3 A 4 (Q 2 )]; this was considered reasonably reliable, because the application of the corresponding extended dPA (90) with M = 2 to the Adler function gave practically the same result for r (D=0) τ [35]. On the other hand, in 2δ AQCD, which was constructed in the (c 2 = −4.9) Lambert renormalization scheme, the two mentioned methods give substantially different results (namely, 0.177 and 0.190, respectively, for the 2δ parameter values given in Table VIII. Since in this scheme the leading-β 0 (LB) parts d  running coupling a(Q 2 ) = α s (Q 2 )/π, have the Padé form where in Eq. (B7c) we denoted In 2δ AQCD, the value of the scheme parameter c 2 = β 2 /β 0 in the Lambert scheme β-function (B7b) can be varied in a restricted interval, −5.9 < c 2 < −2.0 if we require that the pQCD onset scale M 0 is M 0 1 GeV and A max (= A(0)) 1 [24]; the value c 2 = −4.9 is chosen. In 3δ AQCD, the values of the scheme parameters c 2 and c 3 in the LMM β-function (B7c) are adjusted to coincide with the known values of MM scheme [84,85], namely c In Eq. (B9b) the notations ω 1 = c 2 /c 2 1 , ω 2 = c 3 /c 3 1 were used. The squared momentum Q 2 ≡ −q 2 is in general complex, Q 2 = |Q 2 | exp(iφ). The Lambert function W −1 is used when 0 ≤ φ < π, and W +1 when −π ≤ φ < 0. The argument z = z(Q 2 ) in the Lambert functions W ±1 (z) is where Λ L is called the Lambert scale (Λ L 1 GeV 2 ). 21 The solution (B9a) was obtained first in [87], and the solution (B9b) in [86]. The solution (B9a) has one adjustable scheme parameter (c 2 ), and the solution (B9b) has two (c 2 , c 3 ). These solutions, and their spectral functions ρ a (σ) = Ima(Q 2 = −σ − i ), can be evaluated very efficiently with Mathematica 22 , even at high values of |Q 2 | and σ, allowing thus for an efficient evaluation of the dispersion integrals in Eq. (B4). In this work, as the reference point for the underlying coupling constant the value α s (M 2 Z ; MS) = 0.1185 was taken, which then determines the value of the Lambert scale Λ L in the N f = 3 regime of low |Q 2 | 1 GeV 2 . 23 21 In Ref. [86], the expression (B10) without the factor 1/(c 1 e) was used for z, which just redefines the Lambert scale Λ L . 22 In Mathematica, W ∓ (z) is called by the command ProductLog[∓1, z]. 23 The underlying running pQCD coupling a(Q 2 ) in the N f = 3 regime is obtained in the following way: in the four-loop MS scheme the coupling a(M 2 Z ; MS, N f = 5)) = 0.1185/π is RGE-evolved downwards in Q 2 , using the three-loop quark threshold relations [88] at Q 2 = (2mq) 2 (q = b, c; m b = 4.20 GeV and mc = 1.27 GeV), givingā 0 ≡ a(Q 2 0 ; MS) N f =3 = 0.0846346 at Q 2 0 = (2mc) 2 . Then at this scale, the change of the renormalization scheme is made [e.g., cf. Eq. (13) of Ref. [35]], giving a(Q 2 0 ) = 0.0792708 in the c 2 = −4.9 Lambert scheme, and a(Q 2 0 ) = 0.0897919 in the LMM scheme. The five-loop effects in the mentioned RGE evolution are small; namely, in order to reproduce the aforementioned valueā 0 = 0.0846346 when using the polynomial five-loop MS β-function [76] and the corresponding four-loop quark mass thresholds [89], the starting value of αs(M 2 Z ; MS) = 0.118577 is needed which is close to 0.1185. The obtained values of the parameters of the discussed 2δ and 3δ AQCD are given in Table VIII. 24 The scheme parameters c j (j ≥ 2) of the LMM and Lambert schemes, as well as the (five-loop) MS scheme, all at N f = 3, are given in Table IX.
Once a specific coupling A(Q 2 ) is obtained (a → A), the analogs A n (Q 2 ) of the powers a(Q 2 ) n of the underlying pQCD coupling, in general holomorphic AQCD, can be obtained by the construction presented in Ref. [21] for integer n, and in Ref. [39] for general (noninteger) n. The construction of A n (Q 2 ) from A(Q 2 ) for integer n goes via the logarithmic derivatives of A(Q 2 ), in close analogy with the pQCD definitions and relations in Sec. II A.
Here, the construction given in Ref. [21] for integer n will be summarized. Since the coupling A(Q 2 ) is the holomorphic analog of the corresponding underlying pQCD coupling a(Q 2 ) (in the same renormalization scheme), the linearity of the "analytization" a(Q 2 ) → A(Q 2 ) implies that the logarithmic derivatives a n+1 (Q 2 ) of a(Q 2 ), Eq. (2), get replaced (i.e., "analytized") in AQCD by the completely analogous logarithmic derivatives A n+1 (Q 2 ) of A(Q 2 ) where a weak renormalization scale dependence (κ-dependence) now appears due to the truncation effect. The pQCD analog of this expression is the truncated version of the series Eq. (13b), truncated at d N −1 (κ) a N (κQ 2 ). Formally, the truncated series (B12) differs from the full result D(Q 2 ) by a term ∼ A N +1 (∼ a N +1 ∼ a N +1 ). The above truncated series can be rewritten in terms of the coefficients d n (κ) of the original perturbation (power) series (13a) where the power analogs A n are linear combinations of logarithmic derivatives A n+m in complete analogy with the pQCD relations (6) A n = A n + N −n m=1 k n (m) A n+m (n = 2, . . . , N ), where the truncation is performed consistently at A N ; note that A N = A N in this truncation. We recall that the truncated series (B13) has its pQCD analog in the original perturbation (power) series (13a) truncated at d N −1 (κ)a(κQ 2 ) N . 24 Note that the values of the 2δ AQCD parameters differ somewhat from those of Table 2 of Ref. [24] (the case c 2 = −4.9 there), because there the high-momentum reference point was taken to be a(M 2 Z ; MS, N f = 5) = 0.1184/π, while here it is a(M 2 Z ; MS, N f = 5) = 0.1185/π.
We point out that, as long as A(Q 2 ) has some nonperturbative contributions in comparison to its underlying pQCD coupling [such as terms ∼ 1/(Q 2 + M 2 ) k ], we have A n (Q 2 ) = A(Q 2 ) n (n ≥ 2). In such cases, even if the truncation index N in the relations (B14) is very high, we do not have A n (Q 2 ) ≈ A(Q 2 ) n at low values |Q 2 | 1 GeV 2 . 25 In [70] it was argued that if in Eq. (B13) the naive powers A(Q 2 ) n were used instead of A n (Q 2 ), this would bring into the series spurious uncontrollable nonperturbative contributions at |Q 2 | 1 GeV 2 . It is therefore important to use the series in logarithmic derivatives instead, i.e., Eq. (B12) [=Eq. (B13)], and the related extended dPA expressions as explained here in Sec. IV.
In Figs. 7(a),(b) the couplings A(Q 2 ), A 2 (Q 2 ) are presented, as a function of Q 2 > 0, for the considered 2δ and 3δ AQCD, respectively. The corresponding underlying pQCD coupling a(Q 2 ) and the MS couplingā(Q 2 ) are included for comparison (all are for N f = 3). The coupling A 2 (Q 2 ) is obtained by Eq. (B14) with the truncation index N = 4 (i.e., n = 2, N = 4). In these Figures the naive power A(Q 2 ) 2 is included, and we can see clearly that A 2 (Q 2 ) ≈ A(Q 2 ) 2 at low Q 2 . Further, it can be noted that pQCD coupling a(Q 2 ) in the LMM scheme has the branching point at a rather large value Q 2 br = 1.33 GeV 2 , and it is not a pole. In the Lambert scheme, a(Q 2 ) has Q 2 br = 0.068 GeV 2 (it is a pole), and in the MS scheme Q 2 br = 0.42 GeV 2 (it is a pole). These curves can be obtained by using the programs [38] written in Mathematica.
When we are approaching with the e 1 estimate to the true value a 1 , the first term in brackets, (a 1 − e 1 ), becomes smaller, and becomes comparable with the second term at large n (which is ∼ 1/n). This means that we can have, at a large n, cancellation of these two terms. This means that we have a sign change in these differences diff (NLO) (n) at large n when n increases. Therefore, such a sign change at increasing n indicates that we have not sufficiently approximated e 1 to the true value a 1 . On the other hand, when e 1 is adjusted to a 1 with a sufficiently high precision, namely when |a 1 − e 1 | 2|a 2 |/n, the first term (a 1 − e 1 ) in brackets of Eq. (C8) becomes negligible even at large n (we go up to n = 70), and the difference becomes diff (NLO) (n; e 1 ≈ a 1 ) ≈ − S γ 2a 2 This means that we must adjust the value e 1 in such a way that the sequence {diff (NLO) (n; e 1 ); n = 30, . . .} has no sign change (in practice we go up to n = 70). Analogously, the value e 2 must be adjusted so closely to a 2 that the sequence {diff (NNLO) (n; a 1 , e 2 ); n = 30, . . .} has no sign change. These requirements of no sign change represent one of the two mentioned necessary conditions for the determination of a 1 and a 2 , and thus of C 2,k . The other necessary condition is a somewhat arbitrary (although a rather conservative) requirement that at each next order we have at increasing n (40 < n ≤ 70) the velocity of convergence better by at least a factor of two. Specifically, this means that at thus large n, one requires where j = 0, 1, 3, 4. The case when N = j (e.g., when N = 3 = j), the limiting value of the corresponding expression is implied Similarly, the corresponding integrals in the SL case (82b) are (note that 0 < t < 1)